Index: /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp (revision 19739) +++ /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp (revision 19740) @@ -43,5 +43,5 @@ /*Fetch data needed: */ - int hydrology_model; + int hydrology_model,frictionlaw; iomodel->Constant(&hydrology_model,HydrologyModelEnum); @@ -73,5 +73,18 @@ iomodel->FetchDataToInput(elements,HydrologyEnglacialInputEnum); iomodel->FetchDataToInput(elements,HydrologyBumpSpacingEnum); + iomodel->FetchDataToInput(elements,HydrologyBumpHeightEnum); iomodel->FetchDataToInput(elements,HydrologyReynoldsEnum); + iomodel->FetchDataToInput(elements,VxEnum); + iomodel->FetchDataToInput(elements,VyEnum); + + iomodel->Constant(&frictionlaw,FrictionLawEnum); + /*Friction law variables*/ + switch(frictionlaw){ + case 8: + iomodel->FetchDataToInput(elements,FrictionCoefficientEnum); + break; + default: + _error_("Friction law "<< frictionlaw <<" not supported"); + } }/*}}}*/ void HydrologySommersAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ @@ -85,4 +98,5 @@ parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); + parameters->AddObject(iomodel->CopyConstantObject(FrictionLawEnum)); }/*}}}*/ @@ -145,5 +159,7 @@ /*Intermediaries */ IssmDouble Jdet,meltrate,G,dh[2],B,A,n; - IssmDouble gap,bed,thickness,head; + IssmDouble gap,bed,thickness,head,ieb; + IssmDouble lr,br,vx,vy,beta; + IssmDouble alpha2,frictionheat; IssmDouble* xyz_list = NULL; @@ -168,7 +184,15 @@ Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* englacial_input = element->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); + Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); /*Get conductivity from inputs*/ IssmDouble conductivity = GetConductivity(element); + + /*Build friction element, needed later: */ + Friction* friction=new Friction(element,3); /* Start looping on the number of gaussian points: */ @@ -185,4 +209,9 @@ head_input->GetInputValue(&head,gauss); head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss); + englacial_input->GetInputValue(&ieb,gauss); + lr_input->GetInputValue(&lr,gauss); + br_input->GetInputValue(&br,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); /*Get ice A parameter*/ @@ -190,4 +219,16 @@ n_input->GetInputValue(&n,gauss); A=pow(B,-n); + + /*Compute beta term*/ + if(gapGetAlpha2(&alpha2,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + frictionheat=alpha2*(vx*vx+vy*vy); /*Get water and ice pressures*/ @@ -196,5 +237,5 @@ _assert_(pressure_water<=pressure_ice); - meltrate = 1/latentheat*(G+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])); + meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])); _assert_(meltrate>0.); @@ -203,4 +244,6 @@ meltrate*(1/rho_water-1/rho_ice) +A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap + -beta*sqrt(vx*vx+vy*vy) + +ieb )*basis[i]; } @@ -209,4 +252,5 @@ xDelete(xyz_list); xDelete(basis); + delete friction; delete gauss; return pe; @@ -297,2 +341,117 @@ return conductivity; }/*}}}*/ +void HydrologySommersAnalysis::UpdateGapHeight(FemModel* femmodel){/*{{{*/ + + + for(int j=0;jelements->Size();j++){ + Element* element=(Element*)femmodel->elements->GetObjectByOffset(j); + UpdateGapHeight(element); + } + +}/*}}}*/ +void HydrologySommersAnalysis::UpdateGapHeight(Element* element){/*{{{*/ + + /*Skip if water or ice shelf element*/ + if(element->IsFloating()) return; + + /*Intermediaries */ + IssmDouble newgap = 0.; + IssmDouble Jdet,meltrate,G,dh[2],B,A,n,dt; + IssmDouble gap,bed,thickness,head,ieb; + IssmDouble lr,br,vx,vy,beta; + IssmDouble alpha2,frictionheat; + IssmDouble* xyz_list = NULL; + + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&dt,TimesteppingTimeStepEnum); + IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum); + IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input); + Input* head_input = element->GetInput(HydrologyHeadEnum); _assert_(head_input); + Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); + Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* englacial_input = element->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); + Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); + + /*Get conductivity from inputs*/ + IssmDouble conductivity = GetConductivity(element); + + /*Build friction element, needed later: */ + Friction* friction=new Friction(element,3); + + /*Keep track of weights*/ + IssmDouble totalweights=0.; + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + + geothermalflux_input->GetInputValue(&G,gauss); + base_input->GetInputValue(&bed,gauss); + thickness_input->GetInputValue(&thickness,gauss); + gap_input->GetInputValue(&gap,gauss); + head_input->GetInputValue(&head,gauss); + head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss); + englacial_input->GetInputValue(&ieb,gauss); + lr_input->GetInputValue(&lr,gauss); + br_input->GetInputValue(&br,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + + /*Get ice A parameter*/ + B_input->GetInputValue(&B,gauss); + n_input->GetInputValue(&n,gauss); + A=pow(B,-n); + + /*Compute beta term*/ + if(gapGetAlpha2(&alpha2,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + frictionheat=alpha2*(vx*vx+vy*vy); + + /*Get water and ice pressures*/ + IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.); + IssmDouble pressure_water = rho_water*g*(head-bed); + _assert_(pressure_water<=pressure_ice); + + meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])); + _assert_(meltrate>0.); + + newgap += gauss->weight*Jdet*(gap+dt*( + meltrate/rho_ice + -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap + +beta*sqrt(vx*vx+vy*vy) + )); + totalweights +=gauss->weight*Jdet; + } + + /*Divide by connectivity*/ + newgap = newgap/totalweights; + + /*Add new gap as an input*/ + element->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum); + + /*Clean up and return*/ + xDelete(xyz_list); + delete friction; + delete gauss; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.h (revision 19739) +++ /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.h (revision 19740) @@ -33,4 +33,6 @@ /*Intermediaries*/ IssmDouble GetConductivity(Element* element); + void UpdateGapHeight(FemModel* femmodel); + void UpdateGapHeight(Element* element); }; #endif Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 19739) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 19740) @@ -214,6 +214,9 @@ GetAlpha2Coulomb(palpha2,gauss); break; + case 8: + GetAlpha2Sommers(palpha2,gauss); + break; default: - _error_("not supported"); + _error_("Friction law "<< this->law <<" not supported"); } @@ -588,2 +591,35 @@ *palpha2=alpha2; }/*}}}*/ +void Friction::GetAlpha2Sommers(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ + + /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff, with Neff=rho_ice*g*thickness+rho_ice*g*(head-bed)*/ + + /*diverse: */ + IssmDouble pressure_ice,pressure_water; + IssmDouble Neff; + IssmDouble drag_coefficient; + IssmDouble bed,thickness,head; + IssmDouble alpha2; + + /*Recover parameters: */ + element->GetInputValue(&thickness, gauss,ThicknessEnum); + element->GetInputValue(&bed, gauss,BaseEnum); + element->GetInputValue(&head, gauss,HydrologyHeadEnum); + element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + + //From bed and thickness, compute effective pressure when drag is viscous: + pressure_ice = rho_ice*gravity*thickness; + pressure_water = rho_water*gravity*(head-bed); + Neff=pressure_ice-pressure_water; + if(Neff<0.) Neff=0.; + + alpha2=drag_coefficient*drag_coefficient*Neff; + _assert_(!xIsNan(alpha2)); + + /*Assign output pointers:*/ + *palpha2=alpha2; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 19739) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 19740) @@ -40,4 +40,5 @@ void GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2WeertmanTemp(IssmDouble* palpha2,Gauss* gauss); + void GetAlpha2Sommers(IssmDouble* palpha2,Gauss* gauss); }; Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 19739) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 19740) @@ -13,9 +13,8 @@ /*intermediary*/ - int hydrology_model; - bool save_results; - bool modify_loads=true; - bool isefficientlayer; - + int hydrology_model; + bool save_results; + bool modify_loads=true; + bool isefficientlayer; /*first recover parameters common to all solutions*/ @@ -84,4 +83,8 @@ femmodel->SetCurrentConfiguration(HydrologySommersAnalysisEnum); solutionsequence_nonlinear(femmodel,modify_loads); + if(VerboseSolution()) _printf0_(" updating gap height\n"); + HydrologySommersAnalysis* analysis = new HydrologySommersAnalysis(); + analysis->UpdateGapHeight(femmodel); + delete analysis; if(save_results){ Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp (revision 19739) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp (revision 19740) @@ -57,5 +57,5 @@ //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) InputUpdateFromConstantx(femmodel,converged,ConvergedEnum); - InputUpdateFromSolutionx(femmodel,ug); + //InputUpdateFromSolutionx(femmodel,ug); for(;;){