Changeset 20215


Ignore:
Timestamp:
02/19/16 10:25:41 (9 years ago)
Author:
seroussi
Message:

NEW: finished implementing basal melt computation

Location:
issm/trunk-jpl
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r20209 r20215  
    18861886
    18871887        return new GaussPenta(area_coordinates,order_horiz,order_vert);
     1888}
     1889/*}}}*/
     1890Gauss*     Penta::NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){/*{{{*/
     1891        return new GaussPenta(point1,fraction1,fraction2,mainlyfloating,order);
    18881892}
    18891893/*}}}*/
    … …  
    26382642IssmDouble Penta::TotalFloatingBmb(void){/*{{{*/
    26392643
    2640         /*The fbmb[Gt yr-1] of one element is area[m2] * floating basal melt [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
    2641         IssmDouble base,fbmb,rho_ice;
     2644        /*The fbmb[kg yr-1] of one element is area[m2] * melting_rate [kg m^-2 yr^-1]*/
     2645        int        point1;
     2646        bool       mainlyfloating;
     2647        IssmDouble fbmb=0;
     2648        IssmDouble rho_ice,fraction1,fraction2,floatingmelt,Jdet;
    26422649        IssmDouble Total_Fbmb=0;
    26432650        IssmDouble xyz_list[NUMVERTICES][3];
     2651        Gauss*     gauss     = NULL;
     2652
     2653   if(!IsIceInElement() || !IsOnBase())return 0;
    26442654
    26452655        /*Get material parameters :*/
    26462656        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    2647 
    2648         if(!IsIceInElement() || !IsOnBase()) return 0.;
    2649 
     2657        Input* floatingmelt_input = this->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingmelt_input);
     2658        Input* gllevelset_input = this->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
    26502659        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    26512660
    2652         /*First calculate the area of the base (cross section triangle)
    2653          * http://en.wikipedia.org/wiki/Triangle
    2654          * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    2655         base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
    2656 
    2657         /*Now get the average Floating melt rate over the element*/
    2658         Input* fbmb_input = inputs->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fbmb_input);
    2659 
    2660         fbmb_input->GetInputAverage(&fbmb);
    2661         Total_Fbmb=rho_ice*base*fbmb;// floating melt rate on element in kg s-1
     2661        this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
     2662        /* Start  looping on the number of gaussian points: */
     2663        gauss = this->NewGauss(point1,fraction1,fraction2,1-mainlyfloating,3);
     2664        for(int ig=gauss->begin();ig<gauss->end();ig++){
     2665
     2666                gauss->GaussPoint(ig);
     2667                this->JacobianDeterminantBase(&Jdet,&xyz_list[0][0],gauss);
     2668                floatingmelt_input->GetInputValue(&floatingmelt,gauss);
     2669                fbmb+=floatingmelt*Jdet*gauss->weight;
     2670        }
     2671
     2672   Total_Fbmb=rho_ice*fbmb;             // from volume to mass
    26622673
    26632674        /*Return: */
     2675        delete gauss;
    26642676        return Total_Fbmb;
    26652677}
    … …  
    26672679IssmDouble Penta::TotalGroundedBmb(void){/*{{{*/
    26682680
    2669         /*The grounded ice melting rate [Gt yr-1] of one element is area[m2] * grounded melt [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
    2670         IssmDouble base,gbmb,rho_ice;
     2681        /*The gbmb[kg yr-1] of one element is area[m2] * gounded melting rate [kg m^-2 yr^-1]*/
     2682        int        point1;
     2683        bool       mainlyfloating;
     2684        IssmDouble gbmb=0;
     2685        IssmDouble rho_ice,fraction1,fraction2,groundedmelt,Jdet;
    26712686        IssmDouble Total_Gbmb=0;
    26722687        IssmDouble xyz_list[NUMVERTICES][3];
     2688        Gauss*     gauss     = NULL;
     2689
     2690   if(!IsIceInElement() || !IsOnBase())return 0;
    26732691
    26742692        /*Get material parameters :*/
    26752693        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    2676 
    2677         if(!IsIceInElement() || !IsOnBase()) return 0.;
    2678 
     2694        Input* groundedmelt_input = this->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedmelt_input);
     2695        Input* gllevelset_input = this->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
    26792696        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    26802697
    2681         /*First calculate the area of the base (cross section triangle)
    2682          * http://en.wikipedia.org/wiki/Triangle
    2683          * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    2684         base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
    2685 
    2686         /*Now get the average grounded melt over the element*/
    2687         Input* gbmb_input = inputs->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gbmb_input);
    2688 
    2689         gbmb_input->GetInputAverage(&gbmb);
    2690         Total_Gbmb=rho_ice*base*gbmb;// grounded melt on element in kg s-1
     2698        this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
     2699        /* Start  looping on the number of gaussian points: */
     2700        gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,3);
     2701        for(int ig=gauss->begin();ig<gauss->end();ig++){
     2702
     2703                gauss->GaussPoint(ig);
     2704                this->JacobianDeterminantBase(&Jdet,&xyz_list[0][0],gauss);
     2705                groundedmelt_input->GetInputValue(&groundedmelt,gauss);
     2706                gbmb+=groundedmelt*Jdet*gauss->weight;
     2707        }
     2708
     2709   Total_Gbmb=rho_ice*gbmb;             // from volume to mass
    26912710
    26922711        /*Return: */
     2712        delete gauss;
    26932713        return Total_Gbmb;
    26942714}

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r20209 r20215  
    120120                Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
    121121                Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
    122                 Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){_error_("not implemented yet");};
     122                Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order);
    123123                Gauss*         NewGaussBase(int order);
    124124                Gauss*         NewGaussLine(int vertex1,int vertex2,int order);

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r20212 r20215  
    30913091        this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
    30923092        /* Start  looping on the number of gaussian points: */
    3093         gauss = this->NewGauss(point1,fraction1,fraction2,1-mainlyfloating,2);
     3093        gauss = this->NewGauss(point1,fraction1,fraction2,1-mainlyfloating,3);
    30943094        for(int ig=gauss->begin();ig<gauss->end();ig++){
    30953095

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp

    r18521 r20215  
    287287                        }
    288288                        this->weights[ig] = this->weights[ig]*r1*r2;
    289                         this->coords4=xNew<IssmDouble>(numgauss);
    290                         for(ig=0;ig<numgauss;ig++) this->coords4[ig]=-1.0;
    291                 }
     289                }
     290                this->coords4=xNew<IssmDouble>(numgauss);
     291                for(ig=0;ig<numgauss;ig++) this->coords4[ig]=-1.0;
    292292        }
    293293        else{

  • issm/trunk-jpl/test/NightlyRun/test208.py

    r20209 r20215  
    1414md=setmask(md,'all','')
    1515md=parameterize(md,'../Par/SquareShelf.py')
     16md.basalforcings.floatingice_melting_rate[:]=1.
    1617md=setflowequation(md,'SSA','all')
    1718md.cluster=generic('name',oshostname(),'np',3)
    18 md.transient.requested_outputs=['default','GroundedArea','TotalFloatingBmb','TotalGroundedBmb']
     19md.transient.requested_outputs=['default','FloatingArea','GroundedArea','TotalFloatingBmb','TotalGroundedBmb']
    1920md=solve(md,TransientSolutionEnum())
    2021

  • issm/trunk-jpl/test/NightlyRun/test317.m

    r20209 r20215  
    33md=setmask(md,'','');
    44md=parameterize(md,'../Par/SquareSheetConstrained.par');
     5md.basalforcings.groundedice_melting_rate(:)=5;
    56md=extrude(md,3,1.);
    67md=setflowequation(md,'HO','all');

  • issm/trunk-jpl/test/NightlyRun/test317.py

    r20209 r20215  
    1212md=setmask(md,'','')
    1313md=parameterize(md,'../Par/SquareSheetConstrained.py')
     14md.basalforcings.groundedice_melting_rate[:]=5.
    1415md.extrude(3,1.)
    1516md=setflowequation(md,'HO','all')

  • issm/trunk-jpl/test/NightlyRun/test408.m

    r20209 r20215  
    1 %Test Name: SquareSheetShelfTranSSA2d
     1%Test Name: SquareSheetShelfTranSSA3d
    22md=triangle(model(),'../Exp/Square.exp',150000.);
    33md=setmask(md,'../Exp/SquareShelf.exp','');
    44md=parameterize(md,'../Par/SquareSheetShelf.par');
     5md=extrude(md,3,1);
     6md.transient.isthermal=0;
    57md=setflowequation(md,'SSA','all');
    68md.cluster=generic('name',oshostname(),'np',3);
    7 md.transient.requested_outputs={'default','GroundedArea','TotalFloatingBmb','TotalGroundedBmb'};
     9md.transient.requested_outputs={'default','GroundedArea','FloatingArea','TotalFloatingBmb','TotalGroundedBmb','TotalSmb'};
    810md=solve(md,TransientSolutionEnum());
    911
    1012%Fields and tolerances to track changes
    11 field_names     ={'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','GroundedArea1','TotalFloatingBmb1','TotalGroundedBmb1',...
    12         'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','GroundedArea2','TotalFloatingBmb2','TotalGroundedBmb2',...
    13         'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','GroundedArea3','TotalFloatingBmb3','TotalGroundedBmb3'};
    14 field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,...
    15         1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,...
    16         1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13};
     13field_names     ={'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','GroundedArea1','FloatingArea1','TotalFloatingBmb1','TotalGroundedBmb1','TotalSmb1',...
     14        'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','GroundedArea2','FloatingArea2','TotalFloatingBmb2','TotalGroundedBmb2','TotalSmb2',...
     15        'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','GroundedArea3','FloatingArea3','TotalFloatingBmb3','TotalGroundedBmb3','TotalSmb3'};
     16field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,...
     17        1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,...
     18        1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13};
    1719field_values={...
    1820        (md.results.TransientSolution(1).Vx),...
    … …  
    2426        (md.results.TransientSolution(1).Thickness),...
    2527        (md.results.TransientSolution(1).GroundedArea),...
     28        (md.results.TransientSolution(1).FloatingArea),...
    2629        (md.results.TransientSolution(1).TotalFloatingBmb),...
    2730        (md.results.TransientSolution(1).TotalGroundedBmb),...
     31        (md.results.TransientSolution(1).TotalSmb),...
    2832        (md.results.TransientSolution(2).Vx),...
    2933        (md.results.TransientSolution(2).Vy),...
    … …  
    3438        (md.results.TransientSolution(2).Thickness),...
    3539        (md.results.TransientSolution(2).GroundedArea),...
     40        (md.results.TransientSolution(2).FloatingArea),...
    3641        (md.results.TransientSolution(2).TotalFloatingBmb),...
    3742        (md.results.TransientSolution(2).TotalGroundedBmb),...
     43        (md.results.TransientSolution(2).TotalSmb),...
    3844        (md.results.TransientSolution(3).Vx),...
    3945        (md.results.TransientSolution(3).Vy),...
    … …  
    4450        (md.results.TransientSolution(3).Thickness),...
    4551        (md.results.TransientSolution(3).GroundedArea),...
     52        (md.results.TransientSolution(3).FloatingArea),...
    4653        (md.results.TransientSolution(3).TotalFloatingBmb),...
    4754        (md.results.TransientSolution(3).TotalGroundedBmb),...
     55        (md.results.TransientSolution(3).TotalSmb),...
    4856        };

  • issm/trunk-jpl/test/NightlyRun/test408.py

    r20209 r20215  
    1313md=setmask(md,'../Exp/SquareShelf.exp','')
    1414md=parameterize(md,'../Par/SquareSheetShelf.py')
     15md=extrude(md,3,1.)
    1516md=setflowequation(md,'SSA','all')
     17md.transient.isthermal=False
    1618md.cluster=generic('name',oshostname(),'np',3)
    17 md.transient.requested_outputs=['default','GroundedArea','TotalFloatingBmb','TotalGroundedBmb']
     19md.transient.requested_outputs=['default','GroundedArea','FloatingArea','TotalFloatingBmb','TotalGroundedBmb','TotalSmb']
    1820md=solve(md,TransientSolutionEnum())
    1921
    2022#Fields and tolerances to track changes
    21 field_names     =['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','GroundedArea1','TotalFloatingBmb1','TotalGroundedBmb1','Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','GroundedArea2','TotalFloatingBmb2','TotalGroundedBmb2','Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','GroundedArea3','TotalFloatingBmb3','TotalGroundedBmb3']
    22 field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]
     23field_names     =['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','GroundedArea1','TotalFloatingBmb1','TotalGroundedBmb1','TotalSmb1','Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','GroundedArea2','TotalFloatingBmb2','TotalGroundedBmb2','TotalSmb2','Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','GroundedArea3','TotalFloatingBmb3','TotalGroundedBmb3','TotalSmb3']
     24field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]
    2325field_values=[\
    2426        md.results.TransientSolution[0].Vx,\
    … …  
    3032        md.results.TransientSolution[0].Thickness,\
    3133        md.results.TransientSolution[0].GroundedArea,\
     34        md.results.TransientSolution[0].FloatingArea,\
    3235        md.results.TransientSolution[0].TotalFloatingBmb,\
    3336        md.results.TransientSolution[0].TotalGroundedBmb,\
     37        md.results.TransientSolution[0].TotalSmb,\
    3438        md.results.TransientSolution[1].Vx,\
    3539        md.results.TransientSolution[1].Vy,\
    … …  
    4044        md.results.TransientSolution[1].Thickness,\
    4145        md.results.TransientSolution[1].GroundedArea,\
    42         md.results.TransientSolution[2].TotalFloatingBmb,\
    43         md.results.TransientSolution[2].TotalGroundedBmb,\
     46        md.results.TransientSolution[1].FloatingArea,\
     47        md.results.TransientSolution[1].TotalFloatingBmb,\
     48        md.results.TransientSolution[1].TotalGroundedBmb,\
     49        md.results.TransientSolution[1].TotalSmb,\
    4450        md.results.TransientSolution[2].Vx,\
    4551        md.results.TransientSolution[2].Vy,\
    … …  
    5056        md.results.TransientSolution[2].Thickness,\
    5157        md.results.TransientSolution[2].GroundedArea,\
     58        md.results.TransientSolution[2].FloatingArea,\
    5259        md.results.TransientSolution[2].TotalFloatingBmb,\
    5360        md.results.TransientSolution[2].TotalGroundedBmb,\
     61        md.results.TransientSolution[2].TotalSmb,\
    5462        ]

  • issm/trunk-jpl/test/NightlyRun/test433.m

    r19529 r20215  
    4040
    4141md.groundingline.migration='SubelementMigration2';
     42md.transient.isstressbalance=1;
     43md=setflowequation(md,'SSA','all');
    4244md=solve(md,TransientSolutionEnum());
    4345element_on_iceshelf_subelement2=(md.results.TransientSolution.MaskGroundediceLevelset);
Note: See TracChangeset for help on using the changeset viewer.