Changeset 28233

Timestamp:

04/29/24 13:59:08 (11 months ago)

Author:

timhill2

Message:

GlaDS: options to include sheet thickness in hydraulic potential and turn off cavity creep opening

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyGlaDSAnalysis.cpp (modified) (12 diffs)
• classes/Loads/Channel.cpp (modified) (6 diffs)

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

@@ -191 +191 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.omega",HydrologyOmegaEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.istransition",HydrologyIsTransitionEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.isincludesheetthickness",HydrologyIsIncludeSheetThicknessEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.creep_open_flag",HydrologyCreepOpenFlagEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.englacial_void_ratio",HydrologyEnglacialVoidRatioEnum));
 
@@ -241 +243 @@
         /*Get all inputs and parameters*/
         bool istransition;
+        bool isincludesheetthickness;
+        bool creep_open_flag;
         element->FindParam(&istransition,HydrologyIsTransitionEnum);
+        element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
+        element->FindParam(&creep_open_flag,HydrologyCreepOpenFlagEnum);
         IssmDouble alpha     = element->FindParam(HydrologySheetAlphaEnum);
         IssmDouble beta      = element->FindParam(HydrologySheetBetaEnum);
@@ -305 +311 @@
 
                 /*Closing rate term, see Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/
-                phi_0 = rho_water*g*b + rho_ice*g*H;
+                phi_0   = rho_water*g*b + rho_ice*g*H;
+                if(isincludesheetthickness) phi_0 += rho_water*g*h;
                 A=pow(B,-n);
                 v1 = 2./pow(n,n)*A*h*(pow(fabs(phi_0 - phi),n-1.)*( - n));
+                if (!creep_open_flag) {
+                        if (phi_0-phi<0) {
+                                v1 = 0;
+                        }
+                }
                 for(int i=0;i<numnodes;i++){
                         for(int j=0;j<numnodes;j++){
@@ -354 +366 @@
 
         /*Retrieve all inputs and parameters*/
+        bool isincludesheetthickness;
+        bool creep_open_flag;
+        element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
+        element->FindParam(&creep_open_flag,HydrologyCreepOpenFlagEnum);
         element->GetVerticesCoordinates(&xyz_list);
         element->FindParam(&meltflag,HydrologyMeltFlagEnum);
@@ -422 +438 @@
 
                 /*Compute closing rate*/
-                phi_0 = rho_water*g*b + rho_ice*g*H;
+                phi_0   = rho_water*g*b + rho_ice*g*H;
+                if(isincludesheetthickness) phi_0 += rho_water*g*h;
                 A=pow(B,-n);
                 v2 = 2./pow(n,n)*A*h*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
+                if (!creep_open_flag) {
+                        if (phi_0-phi<0) {
+                                v2 = 0.;
+                        }
+                }
 
                 for(int i=0;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(w-v2-m)*basis[i];
@@ -645 +667 @@
 
         /*Retrieve all inputs and parameters*/
+        bool isincludesheetthickness;
+        bool creep_open_flag;
+        element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
+        element->FindParam(&creep_open_flag,HydrologyCreepOpenFlagEnum);
         IssmDouble  dt       = element->FindParam(TimesteppingTimeStepEnum);
         IssmDouble  l_r      = element->FindParam(HydrologyCavitySpacingEnum);
@@ -687 +713 @@
 
                 /*Get values for a few potentials*/
-                phi_0 = rho_water*g*b + rho_ice*g*H;
+                phi_0   = rho_water*g*b + rho_ice*g*H;
+                if(isincludesheetthickness) phi_0 += rho_water*g*h_old;
                 N = phi_0 - phi;
 
@@ -699 +726 @@
                 if(h_old<h_r){
                         alpha = -ub/l_r - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
+                        if (!creep_open_flag) {
+                                if (N<0) {
+                                        alpha = -ub/l_r;
+                                }
+                        }
                         beta  = ub*h_r/l_r;
                 }
                 else{
                         alpha = - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
+                        if (!creep_open_flag) {
+                                if (N<0) {
+                                        alpha = 0.;
+                                }
+                        }
                         beta  = 0.;
                 }
@@ -732 +769 @@
         /*Intermediary*/
         IssmDouble phi_0, phi_m, p_i;
-        IssmDouble H,b,phi;
+        IssmDouble H,b,phi,h;
         IssmDouble oceanLS,iceLS;
 
@@ -738 +775 @@
 
         /*Get thickness and base on nodes to apply cap on water head*/
+        bool isincludesheetthickness;
+        element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
+        Input *h_input       = element->GetInput(HydrologySheetThicknessEnum);    _assert_(h_input);
    IssmDouble* N = xNew<IssmDouble>(numnodes);
         IssmDouble  rho_ice   = element->FindParam(MaterialsRhoIceEnum);
@@ -767 +807 @@
                 oceanLS_input->GetInputValue(&oceanLS,gauss);
                 iceLS_input->GetInputValue(&iceLS,gauss);
+                h_input->GetInputValue(&h,gauss);
 
                 /*Elevation potential*/
@@ -776 +817 @@
                 /*Compute overburden potential*/
                 phi_0 = phi_m + p_i;
+                phi_0   = rho_water*g*b + rho_ice*g*H;
+                if(isincludesheetthickness) phi_0 += rho_water*g*h;
 
                 /*Calculate effective pressure*/

issm/trunk-jpl/src/c/classes/Loads/Channel.cpp

@@ -380 +380 @@
         bool istransition;
         element->FindParam(&istransition,HydrologyIsTransitionEnum);
+        bool isincludesheetthickness;
+        element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
         IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
         IssmDouble mu_water  = element->FindParam(MaterialsMuWaterEnum);
@@ -435 +437 @@
                 /*Get values for a few potentials*/
                 phi_0   = rho_water*g*b + rho_ice*g*H;
+                if(isincludesheetthickness) phi_0 += rho_water*g*h;
                 dphids  = dphi[0]*tx + dphi[1]*ty;
                 dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
@@ -541 +544 @@
         bool istransition;
         element->FindParam(&istransition,HydrologyIsTransitionEnum);
+        bool isincludesheetthickness;
+        element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
         IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
         IssmDouble mu_water  = element->FindParam(MaterialsMuWaterEnum);
@@ -592 +597 @@
                 /*Get values for a few potentials*/
                 phi_0   = rho_water*g*b + rho_ice*g*H;
+                if(isincludesheetthickness) phi_0 += rho_water*g*h;
                 dphids  = dphi[0]*tx + dphi[1]*ty;
                 dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
@@ -691 +697 @@
         bool istransition;
         element->FindParam(&istransition,HydrologyIsTransitionEnum);
+        bool isincludesheetthickness;
+        element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
         IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
@@ -738 +746 @@
         /*Get values for a few potentials*/
         phi_0   = rho_water*g*b + rho_ice*g*H;
+        if(isincludesheetthickness) phi_0 += rho_water*g*h;
         dphids  = dphi[0]*tx + dphi[1]*ty;
         dphimds = rho_water*g*(db[0]*tx + db[1]*ty);