Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 20214) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 20215) @@ -1886,4 +1886,8 @@ return new GaussPenta(area_coordinates,order_horiz,order_vert); +} +/*}}}*/ +Gauss* Penta::NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){/*{{{*/ + return new GaussPenta(point1,fraction1,fraction2,mainlyfloating,order); } /*}}}*/ @@ -2638,28 +2642,36 @@ IssmDouble Penta::TotalFloatingBmb(void){/*{{{*/ - /*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 */ - IssmDouble base,fbmb,rho_ice; + /*The fbmb[kg yr-1] of one element is area[m2] * melting_rate [kg m^-2 yr^-1]*/ + int point1; + bool mainlyfloating; + IssmDouble fbmb=0; + IssmDouble rho_ice,fraction1,fraction2,floatingmelt,Jdet; IssmDouble Total_Fbmb=0; IssmDouble xyz_list[NUMVERTICES][3]; + Gauss* gauss = NULL; + + if(!IsIceInElement() || !IsOnBase())return 0; /*Get material parameters :*/ rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); - - if(!IsIceInElement() || !IsOnBase()) return 0.; - + Input* floatingmelt_input = this->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingmelt_input); + Input* gllevelset_input = this->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input); ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - /*First calculate the area of the base (cross section triangle) - * http://en.wikipedia.org/wiki/Triangle - * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/ - 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])); - - /*Now get the average Floating melt rate over the element*/ - Input* fbmb_input = inputs->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fbmb_input); - - fbmb_input->GetInputAverage(&fbmb); - Total_Fbmb=rho_ice*base*fbmb;// floating melt rate on element in kg s-1 + this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + /* Start looping on the number of gaussian points: */ + gauss = this->NewGauss(point1,fraction1,fraction2,1-mainlyfloating,3); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + this->JacobianDeterminantBase(&Jdet,&xyz_list[0][0],gauss); + floatingmelt_input->GetInputValue(&floatingmelt,gauss); + fbmb+=floatingmelt*Jdet*gauss->weight; + } + + Total_Fbmb=rho_ice*fbmb; // from volume to mass /*Return: */ + delete gauss; return Total_Fbmb; } @@ -2667,28 +2679,36 @@ IssmDouble Penta::TotalGroundedBmb(void){/*{{{*/ - /*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 */ - IssmDouble base,gbmb,rho_ice; + /*The gbmb[kg yr-1] of one element is area[m2] * gounded melting rate [kg m^-2 yr^-1]*/ + int point1; + bool mainlyfloating; + IssmDouble gbmb=0; + IssmDouble rho_ice,fraction1,fraction2,groundedmelt,Jdet; IssmDouble Total_Gbmb=0; IssmDouble xyz_list[NUMVERTICES][3]; + Gauss* gauss = NULL; + + if(!IsIceInElement() || !IsOnBase())return 0; /*Get material parameters :*/ rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); - - if(!IsIceInElement() || !IsOnBase()) return 0.; - + Input* groundedmelt_input = this->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedmelt_input); + Input* gllevelset_input = this->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input); ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - /*First calculate the area of the base (cross section triangle) - * http://en.wikipedia.org/wiki/Triangle - * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/ - 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])); - - /*Now get the average grounded melt over the element*/ - Input* gbmb_input = inputs->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gbmb_input); - - gbmb_input->GetInputAverage(&gbmb); - Total_Gbmb=rho_ice*base*gbmb;// grounded melt on element in kg s-1 + this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + /* Start looping on the number of gaussian points: */ + gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,3); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + this->JacobianDeterminantBase(&Jdet,&xyz_list[0][0],gauss); + groundedmelt_input->GetInputValue(&groundedmelt,gauss); + gbmb+=groundedmelt*Jdet*gauss->weight; + } + + Total_Gbmb=rho_ice*gbmb; // from volume to mass /*Return: */ + delete gauss; return Total_Gbmb; } Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 20214) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 20215) @@ -120,5 +120,5 @@ Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");}; Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert); - Gauss* NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){_error_("not implemented yet");}; + Gauss* NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order); Gauss* NewGaussBase(int order); Gauss* NewGaussLine(int vertex1,int vertex2,int order); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20214) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20215) @@ -3091,5 +3091,5 @@ this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); /* Start looping on the number of gaussian points: */ - gauss = this->NewGauss(point1,fraction1,fraction2,1-mainlyfloating,2); + gauss = this->NewGauss(point1,fraction1,fraction2,1-mainlyfloating,3); for(int ig=gauss->begin();igend();ig++){ Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 20214) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 20215) @@ -287,7 +287,7 @@ } this->weights[ig] = this->weights[ig]*r1*r2; - this->coords4=xNew(numgauss); - for(ig=0;igcoords4[ig]=-1.0; - } + } + this->coords4=xNew(numgauss); + for(ig=0;igcoords4[ig]=-1.0; } else{ Index: /issm/trunk-jpl/test/NightlyRun/test208.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test208.py (revision 20214) +++ /issm/trunk-jpl/test/NightlyRun/test208.py (revision 20215) @@ -14,7 +14,8 @@ md=setmask(md,'all','') md=parameterize(md,'../Par/SquareShelf.py') +md.basalforcings.floatingice_melting_rate[:]=1. md=setflowequation(md,'SSA','all') md.cluster=generic('name',oshostname(),'np',3) -md.transient.requested_outputs=['default','GroundedArea','TotalFloatingBmb','TotalGroundedBmb'] +md.transient.requested_outputs=['default','FloatingArea','GroundedArea','TotalFloatingBmb','TotalGroundedBmb'] md=solve(md,TransientSolutionEnum()) Index: /issm/trunk-jpl/test/NightlyRun/test317.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test317.m (revision 20214) +++ /issm/trunk-jpl/test/NightlyRun/test317.m (revision 20215) @@ -3,4 +3,5 @@ md=setmask(md,'',''); md=parameterize(md,'../Par/SquareSheetConstrained.par'); +md.basalforcings.groundedice_melting_rate(:)=5; md=extrude(md,3,1.); md=setflowequation(md,'HO','all'); Index: /issm/trunk-jpl/test/NightlyRun/test317.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test317.py (revision 20214) +++ /issm/trunk-jpl/test/NightlyRun/test317.py (revision 20215) @@ -12,4 +12,5 @@ md=setmask(md,'','') md=parameterize(md,'../Par/SquareSheetConstrained.py') +md.basalforcings.groundedice_melting_rate[:]=5. md.extrude(3,1.) md=setflowequation(md,'HO','all') Index: /issm/trunk-jpl/test/NightlyRun/test408.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test408.m (revision 20214) +++ /issm/trunk-jpl/test/NightlyRun/test408.m (revision 20215) @@ -1,18 +1,20 @@ -%Test Name: SquareSheetShelfTranSSA2d +%Test Name: SquareSheetShelfTranSSA3d md=triangle(model(),'../Exp/Square.exp',150000.); md=setmask(md,'../Exp/SquareShelf.exp',''); md=parameterize(md,'../Par/SquareSheetShelf.par'); +md=extrude(md,3,1); +md.transient.isthermal=0; md=setflowequation(md,'SSA','all'); md.cluster=generic('name',oshostname(),'np',3); -md.transient.requested_outputs={'default','GroundedArea','TotalFloatingBmb','TotalGroundedBmb'}; +md.transient.requested_outputs={'default','GroundedArea','FloatingArea','TotalFloatingBmb','TotalGroundedBmb','TotalSmb'}; md=solve(md,TransientSolutionEnum()); %Fields and tolerances to track changes -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'}; -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}; +field_names ={'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','GroundedArea1','FloatingArea1','TotalFloatingBmb1','TotalGroundedBmb1','TotalSmb1',... + 'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','GroundedArea2','FloatingArea2','TotalFloatingBmb2','TotalGroundedBmb2','TotalSmb2',... + 'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','GroundedArea3','FloatingArea3','TotalFloatingBmb3','TotalGroundedBmb3','TotalSmb3'}; +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,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; field_values={... (md.results.TransientSolution(1).Vx),... @@ -24,6 +26,8 @@ (md.results.TransientSolution(1).Thickness),... (md.results.TransientSolution(1).GroundedArea),... + (md.results.TransientSolution(1).FloatingArea),... (md.results.TransientSolution(1).TotalFloatingBmb),... (md.results.TransientSolution(1).TotalGroundedBmb),... + (md.results.TransientSolution(1).TotalSmb),... (md.results.TransientSolution(2).Vx),... (md.results.TransientSolution(2).Vy),... @@ -34,6 +38,8 @@ (md.results.TransientSolution(2).Thickness),... (md.results.TransientSolution(2).GroundedArea),... + (md.results.TransientSolution(2).FloatingArea),... (md.results.TransientSolution(2).TotalFloatingBmb),... (md.results.TransientSolution(2).TotalGroundedBmb),... + (md.results.TransientSolution(2).TotalSmb),... (md.results.TransientSolution(3).Vx),... (md.results.TransientSolution(3).Vy),... @@ -44,5 +50,7 @@ (md.results.TransientSolution(3).Thickness),... (md.results.TransientSolution(3).GroundedArea),... + (md.results.TransientSolution(3).FloatingArea),... (md.results.TransientSolution(3).TotalFloatingBmb),... (md.results.TransientSolution(3).TotalGroundedBmb),... + (md.results.TransientSolution(3).TotalSmb),... }; Index: /issm/trunk-jpl/test/NightlyRun/test408.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test408.py (revision 20214) +++ /issm/trunk-jpl/test/NightlyRun/test408.py (revision 20215) @@ -13,12 +13,14 @@ md=setmask(md,'../Exp/SquareShelf.exp','') md=parameterize(md,'../Par/SquareSheetShelf.py') +md=extrude(md,3,1.) md=setflowequation(md,'SSA','all') +md.transient.isthermal=False md.cluster=generic('name',oshostname(),'np',3) -md.transient.requested_outputs=['default','GroundedArea','TotalFloatingBmb','TotalGroundedBmb'] +md.transient.requested_outputs=['default','GroundedArea','FloatingArea','TotalFloatingBmb','TotalGroundedBmb','TotalSmb'] md=solve(md,TransientSolutionEnum()) #Fields and tolerances to track changes -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'] -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] +field_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'] +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,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13] field_values=[\ md.results.TransientSolution[0].Vx,\ @@ -30,6 +32,8 @@ md.results.TransientSolution[0].Thickness,\ md.results.TransientSolution[0].GroundedArea,\ + md.results.TransientSolution[0].FloatingArea,\ md.results.TransientSolution[0].TotalFloatingBmb,\ md.results.TransientSolution[0].TotalGroundedBmb,\ + md.results.TransientSolution[0].TotalSmb,\ md.results.TransientSolution[1].Vx,\ md.results.TransientSolution[1].Vy,\ @@ -40,6 +44,8 @@ md.results.TransientSolution[1].Thickness,\ md.results.TransientSolution[1].GroundedArea,\ - md.results.TransientSolution[2].TotalFloatingBmb,\ - md.results.TransientSolution[2].TotalGroundedBmb,\ + md.results.TransientSolution[1].FloatingArea,\ + md.results.TransientSolution[1].TotalFloatingBmb,\ + md.results.TransientSolution[1].TotalGroundedBmb,\ + md.results.TransientSolution[1].TotalSmb,\ md.results.TransientSolution[2].Vx,\ md.results.TransientSolution[2].Vy,\ @@ -50,5 +56,7 @@ md.results.TransientSolution[2].Thickness,\ md.results.TransientSolution[2].GroundedArea,\ + md.results.TransientSolution[2].FloatingArea,\ md.results.TransientSolution[2].TotalFloatingBmb,\ md.results.TransientSolution[2].TotalGroundedBmb,\ + md.results.TransientSolution[2].TotalSmb,\ ] Index: /issm/trunk-jpl/test/NightlyRun/test433.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test433.m (revision 20214) +++ /issm/trunk-jpl/test/NightlyRun/test433.m (revision 20215) @@ -40,4 +40,6 @@ md.groundingline.migration='SubelementMigration2'; +md.transient.isstressbalance=1; +md=setflowequation(md,'SSA','all'); md=solve(md,TransientSolutionEnum()); element_on_iceshelf_subelement2=(md.results.TransientSolution.MaskGroundediceLevelset);