Changeset 26137

Timestamp:

03/23/21 13:56:17 (4 years ago)

Author:

tsantos

Message:

CHG: working on mono layer HO

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (modified) (5 diffs)

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

@@ -2555 +2555 @@
 
         /*Intermediaries*/
-        IssmDouble  Jdet,thickness,bed,sealevel,water_pressure,ice_pressure;
+        IssmDouble  Jdet,thickness,base,sealevel,water_pressure,ice_pressure;
         IssmDouble  surface_under_water,base_under_water,pressure;
         IssmDouble *xyz_list = NULL;
@@ -2583 +2583 @@
         while(gauss->next()){
                 thickness_input->GetInputValue(&thickness,gauss);
-                base_input->GetInputValue(&bed,gauss);
+                base_input->GetInputValue(&base,gauss);
                 sealevel_input->GetInputValue(&sealevel,gauss);
                 element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
                 element->NodalFunctions(basis,gauss);
 
-                surface_under_water = min(0.,thickness+bed-sealevel); // 0 if the top of the glacier is above water level
-                base_under_water    = min(0.,bed-sealevel);           // 0 if the bottom of the glacier is above water level
+                surface_under_water = min(0.,thickness+base-sealevel); // 0 if the top of the glacier is above water level
+                base_under_water    = min(0.,base-sealevel);           // 0 if the bottom of the glacier is above water level
                 water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
                 ice_pressure   = 1.0/2.0*gravity*rho_ice*thickness*thickness;
@@ -2949 +2949 @@
 
         /*Intermediaries*/
-        IssmDouble  Jdet,thickness,bed,sealevel,water_pressure,ice_pressure;
+        IssmDouble  Jdet,thickness,base,sealevel,water_pressure,ice_pressure;
+        IssmDouble  water_pressure_sh,ice_pressure_sh,pressure_sh,n,s,b;
         IssmDouble  surface_under_water,base_under_water,pressure;
         IssmDouble *xyz_list = NULL;
@@ -2964 +2965 @@
         /*Retrieve all inputs and parameters*/
         Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
-        Input* base_input       = element->GetInput(BaseEnum);       _assert_(base_input);
-        Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
+        Input* base_input      = element->GetInput(BaseEnum);       _assert_(base_input);
+        Input* sealevel_input  = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
+        Input* n_input         = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
         IssmDouble rho_water   = element->FindParam(MaterialsRhoSeawaterEnum);
         IssmDouble rho_ice     = element->FindParam(MaterialsRhoIceEnum);
@@ -2977 +2979 @@
         while(gauss->next()){
                 thickness_input->GetInputValue(&thickness,gauss);
-                base_input->GetInputValue(&bed,gauss);
+                base_input->GetInputValue(&base,gauss);
+                n_input->GetInputValue(&n,gauss);
                 sealevel_input->GetInputValue(&sealevel,gauss);
                 element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
                 element->NodalFunctions(basis,gauss);
 
-                surface_under_water = min(0.,thickness+bed-sealevel); // 0 if the top of the glacier is above water level
-                base_under_water    = min(0.,bed-sealevel);           // 0 if the bottom of the glacier is above water level
+                surface_under_water = min(0.,thickness+base-sealevel); // 0 if the top of the glacier is above water level
+                base_under_water    = min(0.,base-sealevel);           // 0 if the bottom of the glacier is above water level
+                /*Vertically integrated pressure - SSA type*/
                 water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
                 ice_pressure   = 1.0/2.0*gravity*rho_ice*thickness*thickness;
                 pressure = ice_pressure + water_pressure;
+                /*Vertically integrated pressure - HO type*/
+                s=max(0.,thickness+base-sealevel);
+                b=min(0.,base-sealevel); 
+                water_pressure_sh = gravity*rho_water*(-b*b/2 + pow(-b,n+2)*( -s/(n+2) -b/(n+3) )/pow(thickness,n+1));
+                ice_pressure_sh   = gravity*rho_ice*thickness*thickness*(n+1)/(2*(n+3));
+                pressure_sh = ice_pressure_sh + water_pressure_sh;
 
                 for (int i=0;i<numnodes;i++){
-                        pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
-                        pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
+                        pe->values[i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; // F1
+                        pe->values[i+3]+= pressure_sh*Jdet*gauss->weight*normal[0]*basis[i]; // F2
+                        pe->values[i+6]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; // F3
+                        pe->values[i+9]+= pressure_sh*Jdet*gauss->weight*normal[1]*basis[i]; // F4
                 }
         }