Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23940) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23941) @@ -113,9 +113,11 @@ ElementMatrix* HydrologyGlaDSAnalysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("Not implemented"); - /*Intermediaries */ - IssmDouble Jdet; + IssmDouble Jdet,dphi[3],h,k; IssmDouble* xyz_list = NULL; + + /*Hard coded coefficients*/ + const IssmPDouble alpha = 5./4.; + const IssmPDouble beta = 3./2.; /*Fetch number of nodes and dof for this finite element*/ @@ -130,11 +132,14 @@ element->GetVerticesCoordinates(&xyz_list); - /*Get englacial storage coefficient*/ - IssmDouble storage,dt; - element->FindParam(&storage,HydrologyStorageEnum); + /*Get all inputs and parameters*/ + IssmDouble dt,c_t; element->FindParam(&dt,TimesteppingTimeStepEnum); + element->FindParam(&c_t,HydrologyPressureMeltCoefficientEnum); + element->FindParam(&k,HydrologySheetConductivityEnum); + Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); + Input* h_input = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input); /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(1); + Gauss* gauss=element->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); @@ -144,8 +149,20 @@ element->NodalFunctions(basis,gauss); + phi_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss); + h_input->GetInputValue(&h,gauss); + + /*Get norm of gradient of hydraulic potential and make sure it is >0*/ + IssmDouble normgradphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]); + if(normgradphi < 1.e-12) normgradphi = 1.e-12; + + IssmDouble coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.); + for(int i=0;ivalues[i*numnodes+j] += gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]) - + gauss->weight*Jdet*storage/dt*basis[i]*basis[j]; + + /*Diffusive term*/ + Ke->values[i*numnodes+j] += gauss->weight*Jdet*( + coeff*dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + + coeff*dbasis[1*numnodes+i]*dbasis[1*numnodes+j]); } } @@ -164,12 +181,8 @@ if(element->IsFloating()) return NULL; - _error_("not implemented"); - /*Intermediaries */ - IssmDouble Jdet,meltrate,G,dh[2],B,A,n; - IssmDouble gap,bed,thickness,head,ieb,head_old; - IssmDouble lr,br,vx,vy,beta; - IssmDouble alpha2,frictionheat; - IssmDouble PMPheat,dpressure_water[2],dbed[2]; + IssmDouble Jdet,w,v,vx,vy,ub,h,N; + IssmDouble G,m,frictionheat,alpha2; + IssmDouble A,B,n; IssmDouble* xyz_list = NULL; @@ -183,27 +196,15 @@ /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); - IssmDouble latentheat = element->FindParam(MaterialsLatentheatEnum); - IssmDouble g = element->FindParam(ConstantsGEnum); - IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); - IssmDouble rho_water = element->FindParam(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); - Input* headold_input = element->GetInput(HydrologyHeadOldEnum); _assert_(headold_input); - - - /*Get englacial storage coefficient*/ - IssmDouble storage,dt; - element->FindParam(&storage,HydrologyStorageEnum); - element->FindParam(&dt,TimesteppingTimeStepEnum); + IssmDouble h_r = element->FindParam(HydrologyBumpHeightEnum); + IssmDouble l_r = element->FindParam(HydrologyBumpSpacingEnum); + IssmDouble L = element->FindParam(MaterialsLatentheatEnum); + IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); + Input* vx_input = element->GetInput(VxEnum);_assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum);_assert_(vy_input); + Input* N_input = element->GetInput(EffectivePressureEnum); _assert_(N_input); + Input* h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input); + Input* G_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(G_input); + Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); /*Build friction element, needed later: */ @@ -217,71 +218,46 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis,gauss); - geothermalflux_input->GetInputValue(&G,gauss); - base_input->GetInputValue(&bed,gauss); - base_input->GetInputDerivativeValue(&dbed[0],xyz_list,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); + + /*Get input values at gauss points*/ vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); - headold_input->GetInputValue(&head_old,gauss); - - /*Get ice A parameter*/ + h_input->GetInputValue(&h,gauss); + G_input->GetInputValue(&G,gauss); B_input->GetInputValue(&B,gauss); n_input->GetInputValue(&n,gauss); - A=pow(B,-n); - - /*Compute beta term*/ - if(gapGetInputValue(&N,gauss); + + /*Get basal velocity*/ + ub = sqrt(vx*vx + vy*vy); + + /*Compute cavity opening w*/ + w = 0.; + if(hGetAlpha2(&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); - if(pressure_water>pressure_ice) pressure_water = pressure_ice; - - /*Get water pressure from previous time step to use in lagged creep term*/ - IssmDouble pressure_water_old = rho_water*g*(head_old-bed); - if(pressure_water_old>pressure_ice) pressure_water_old = pressure_ice; - - /*Compute change in sensible heat due to changes in pressure melting point*/ - dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]); - dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]); - PMPheat=0.; - - meltrate = 1/latentheat*(G+frictionheat+rho_water*g*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat); - _assert_(meltrate>0.); - - for(int i=0;ivalues[i]+=Jdet*gauss->weight* - ( - 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 - +storage*head_old/dt - )*basis[i]; - } + frictionheat=alpha2*ub*ub; + + /*Compute melt*/ + m = (G + frictionheat)/(rho_ice*L); + + /*Compute closing rate*/ + A=pow(B,-n); + v = 2./pow(n,n)*A*h*pow(fabs(N),n-1.)*N; + + for(int i=0;ivalues[i]+= - Jdet*gauss->weight*(w-v-m)*basis[i]; + } + /*Clean up and return*/ xDelete(xyz_list); xDelete(basis); + delete gauss; delete friction; - delete gauss; return pe; }/*}}}*/ void HydrologyGlaDSAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ - _error_("not implemented"); - element->GetSolutionFromInputsOneDof(solution,HydrologyHeadEnum); + + element->GetSolutionFromInputsOneDof(solution,HydraulicPotentialEnum); + }/*}}}*/ void HydrologyGlaDSAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ @@ -290,13 +266,7 @@ void HydrologyGlaDSAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - _error_("not implemented"); - /*Intermediary*/ - IssmDouble dh[3]; int* doflist = NULL; - IssmDouble* xyz_list = NULL; - - /*Get gravity from parameters*/ - IssmDouble g = element->FindParam(ConstantsGEnum); + IssmDouble phi_0, phi_m, pi; /*Fetch number of nodes for this finite element*/ @@ -308,17 +278,13 @@ /*Get thickness and base on nodes to apply cap on water head*/ - IssmDouble* eff_pressure = xNew(numnodes); + IssmDouble* N = xNew(numnodes); + IssmDouble* h = xNew(numnodes); IssmDouble* thickness = xNew(numnodes); IssmDouble* bed = xNew(numnodes); IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); + IssmDouble g = element->FindParam(ConstantsGEnum); element->GetInputListOnNodes(&thickness[0],ThicknessEnum); element->GetInputListOnNodes(&bed[0],BaseEnum); - - /*Get head from previous time-step and under-relaxation coefficient to use in under-relaxation for nonlinear convergence*/ - IssmDouble* head_old = xNew(numnodes); - element->GetInputListOnNodes(&head_old[0],HydrologyHeadEnum); - IssmDouble relaxation; - element->FindParam(&relaxation,HydrologyRelaxationEnum); /*Use the dof list to index into the solution vector: */ @@ -326,19 +292,15 @@ values[i]=solution[doflist[i]]; - /*make sure that p_water hrho_ice*thickness[i]/rho_water+bed[i]){ - values[i] = rho_ice*thickness[i]/rho_water+bed[i]; - } - - /*Make sure that negative pressure is not allowed*/ - // if(values[i](values[i])) _error_("NaN found in solution vector"); @@ -347,20 +309,14 @@ /*Add input to the element: */ - element->AddInput(HydrologyHeadEnum,values,element->GetElementType()); - element->AddInput(EffectivePressureEnum,eff_pressure,P1Enum); - - /*Update reynolds number according to new solution*/ - element->GetVerticesCoordinates(&xyz_list); - Input* head_input = element->GetInput(HydrologyHeadEnum);_assert_(head_input); - head_input->GetInputDerivativeAverageValue(&dh[0],xyz_list); + element->AddInput(HydraulicPotentialEnum,values,element->GetElementType()); + element->AddInput(EffectivePressureEnum,N,P1Enum); /*Free resources:*/ + xDelete(doflist); xDelete(values); + xDelete(N); + xDelete(h); + xDelete(bed); xDelete(thickness); - xDelete(bed); - xDelete(xyz_list); - xDelete(doflist); - xDelete(eff_pressure); - xDelete(head_old); }/*}}}*/ void HydrologyGlaDSAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 23940) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 23941) @@ -355,5 +355,5 @@ /*Get gravity from parameters*/ - IssmDouble g = element->FindParam(ConstantsGEnum); + IssmDouble g = element->FindParam(ConstantsGEnum); /*Fetch number of nodes for this finite element*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23940) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23941) @@ -536,5 +536,7 @@ GradientEnum, GroundinglineHeightEnum, + HydraulicPotentialEnum, HydrologyBasalFluxEnum, + HydrologySheetThicknessEnum, HydrologyBumpHeightEnum, HydrologyBumpSpacingEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23940) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23941) @@ -542,5 +542,7 @@ case GradientEnum : return "Gradient"; case GroundinglineHeightEnum : return "GroundinglineHeight"; + case HydraulicPotentialEnum : return "HydraulicPotential"; case HydrologyBasalFluxEnum : return "HydrologyBasalFlux"; + case HydrologySheetThicknessEnum : return "HydrologySheetThickness"; case HydrologyBumpHeightEnum : return "HydrologyBumpHeight"; case HydrologyBumpSpacingEnum : return "HydrologyBumpSpacing"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23940) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23941) @@ -554,5 +554,7 @@ else if (strcmp(name,"Gradient")==0) return GradientEnum; else if (strcmp(name,"GroundinglineHeight")==0) return GroundinglineHeightEnum; + else if (strcmp(name,"HydraulicPotential")==0) return HydraulicPotentialEnum; else if (strcmp(name,"HydrologyBasalFlux")==0) return HydrologyBasalFluxEnum; + else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum; else if (strcmp(name,"HydrologyBumpHeight")==0) return HydrologyBumpHeightEnum; else if (strcmp(name,"HydrologyBumpSpacing")==0) return HydrologyBumpSpacingEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; else if (strcmp(name,"RheologyBInitialguessMisfit")==0) return RheologyBInitialguessMisfitEnum; - else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum; - else if (strcmp(name,"Sealevel")==0) return SealevelEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"SealevelEustaticMask")==0) return SealevelEustaticMaskEnum; + if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum; + else if (strcmp(name,"Sealevel")==0) return SealevelEnum; + else if (strcmp(name,"SealevelEustaticMask")==0) return SealevelEustaticMaskEnum; else if (strcmp(name,"SealevelriseCumDeltathickness")==0) return SealevelriseCumDeltathicknessEnum; else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum; @@ -750,10 +752,10 @@ else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum; else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; - else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum; - else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; + if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum; + else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum; + else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum; else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum; @@ -873,10 +875,10 @@ else if (strcmp(name,"Outputdefinition67")==0) return Outputdefinition67Enum; else if (strcmp(name,"Outputdefinition68")==0) return Outputdefinition68Enum; - else if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum; - else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum; + if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum; + else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum; + else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum; else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum; else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum; @@ -996,10 +998,10 @@ else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum; else if (strcmp(name,"FemModel")==0) return FemModelEnum; - else if (strcmp(name,"FileParam")==0) return FileParamEnum; - else if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum; + if (strcmp(name,"FileParam")==0) return FileParamEnum; + else if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum; + else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum; else if (strcmp(name,"FloatingAreaScaled")==0) return FloatingAreaScaledEnum; else if (strcmp(name,"FloatingMeltRate")==0) return FloatingMeltRateEnum; @@ -1119,10 +1121,10 @@ else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum; else if (strcmp(name,"MINI")==0) return MINIEnum; - else if (strcmp(name,"MinVel")==0) return MinVelEnum; - else if (strcmp(name,"MinVx")==0) return MinVxEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"MinVy")==0) return MinVyEnum; + if (strcmp(name,"MinVel")==0) return MinVelEnum; + else if (strcmp(name,"MinVx")==0) return MinVxEnum; + else if (strcmp(name,"MinVy")==0) return MinVyEnum; else if (strcmp(name,"MinVz")==0) return MinVzEnum; else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum; @@ -1242,10 +1244,10 @@ else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; - else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; - else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"Tria")==0) return TriaEnum; + if (strcmp(name,"TransientParam")==0) return TransientParamEnum; + else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; + else if (strcmp(name,"Tria")==0) return TriaEnum; else if (strcmp(name,"TriaInput")==0) return TriaInputEnum; else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum;