Changeset 23941

Timestamp:

05/29/19 10:14:47 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: worked on RHS and K of GlaDS

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyGlaDSAnalysis.cpp (modified) (10 diffs)
• analyses/HydrologyShaktiAnalysis.cpp (modified) (1 diff)
• shared/Enum/EnumDefinitions.h (modified) (1 diff)
• shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• shared/Enum/StringToEnumx.cpp (modified) (7 diffs)

issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp ¶

@@ -113 +113 @@
 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 +132 @@
         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();ig<gauss->end();ig++){
                 gauss->GaussPoint(ig);
@@ -144 +149 @@
                 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;i<numnodes;i++){
                         for(int j=0;j<numnodes;j++){
-                                Ke->values[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 +181 @@
         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 +196 @@
         /*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 +218 @@
                 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(gap<br)
-                 beta = (br-gap)/lr;
-                else
-                 beta = 0.;
+                N_input->GetInputValue(&N,gauss);
+
+                /*Get basal velocity*/
+                ub = sqrt(vx*vx + vy*vy);
+
+                /*Compute cavity opening w*/
+                w  = 0.;
+                if(h<h_r) w = ub*(h_r-h)/l_r;
 
                 /*Compute frictional heat flux*/
                 friction->GetAlpha2(&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;i<numnodes;i++) pe->values[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;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(w-v-m)*basis[i];
+        }
+
         /*Clean up and return*/
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(basis);
+        delete gauss;
         delete friction;
-        delete gauss;
         return pe;
 }/*}}}*/
 void           HydrologyGlaDSAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-        _error_("not implemented");
-        element->GetSolutionFromInputsOneDof(solution,HydrologyHeadEnum);
+
+        element->GetSolutionFromInputsOneDof(solution,HydraulicPotentialEnum);
+
 }/*}}}*/
 void           HydrologyGlaDSAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
@@ -290 +266 @@
 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 +278 @@
 
         /*Get thickness and base on nodes to apply cap on water head*/
-   IssmDouble* eff_pressure = xNew<IssmDouble>(numnodes);
+   IssmDouble* N = xNew<IssmDouble>(numnodes);
+        IssmDouble* h = xNew<IssmDouble>(numnodes);
         IssmDouble* thickness = xNew<IssmDouble>(numnodes);
         IssmDouble* bed       = xNew<IssmDouble>(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<IssmDouble>(numnodes); 
-        element->GetInputListOnNodes(&head_old[0],HydrologyHeadEnum);
-   IssmDouble relaxation; 
-        element->FindParam(&relaxation,HydrologyRelaxationEnum);
 
         /*Use the dof list to index into the solution vector: */
@@ -326 +292 @@
                 values[i]=solution[doflist[i]];
 
-                /*make sure that p_water<p_ice ->  h<rho_i H/rho_w + zb*/
-                if(values[i]>rho_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]<bed[i]){
-        //              values[i] = bed[i];
-        //      }
-
-                /*Under-relaxation*/
-           values[i] = head_old[i] - relaxation*(head_old[i]-values[i]);
+                /*Elevation potential*/
+                phi_m = rho_water*g*bed[i];
+
+                /*Compute overburden pressure*/
+                pi = rho_ice*g*thickness[i];
+
+                /*Copmute overburden potential*/
+                phi_0 = phi_m + pi;
 
                 /*Calculate effective pressure*/
-                eff_pressure[i] = rho_ice*g*thickness[i] - rho_water*g*(values[i]-bed[i]);
+                N[i] = phi_0 - values[i];
 
                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
@@ -347 +309 @@
 
         /*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<int>(doflist);
         xDelete<IssmDouble>(values);
+        xDelete<IssmDouble>(N);
+        xDelete<IssmDouble>(h);
+        xDelete<IssmDouble>(bed);
         xDelete<IssmDouble>(thickness);
-        xDelete<IssmDouble>(bed);
-        xDelete<IssmDouble>(xyz_list);
-        xDelete<int>(doflist);
-        xDelete<IssmDouble>(eff_pressure);
-   xDelete<IssmDouble>(head_old);
 }/*}}}*/
 void           HydrologyGlaDSAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/

issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp ¶

@@ -355 +355 @@
 
         /*Get gravity from parameters*/
-           IssmDouble  g = element->FindParam(ConstantsGEnum);
+        IssmDouble  g = element->FindParam(ConstantsGEnum);
 
         /*Fetch number of nodes for this finite element*/

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h ¶

@@ -536 +536 @@
         GradientEnum,
         GroundinglineHeightEnum,
+        HydraulicPotentialEnum,
         HydrologyBasalFluxEnum,
+        HydrologySheetThicknessEnum,
         HydrologyBumpHeightEnum,
         HydrologyBumpSpacingEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp ¶

@@ -542 +542 @@
                 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";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp ¶

@@ -554 +554 @@
               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 +629 @@
               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 +752 @@
               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 +875 @@
               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 +998 @@
               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 +1121 @@
               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 +1244 @@
               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;