Changeset 21145

Timestamp:

08/17/16 09:52:54 (9 years ago)

Author:

aleahsommers

Message:

CHG: Added effective pressure as an output, added term in meltrate to include changes in pressure melting point

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologySommersAnalysis.cpp (modified) (17 diffs)
• cores/hydrology_core.cpp (modified) (1 diff)

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

@@ -8 +8 @@
 #define OMEGA 0.001    // parameter controlling transition to nonlinear resistance in basal system (dimensionless)
 #define NU    1.787e-6 //kinematic water viscosity m^2/s
+#define CT    7.5e-8  // Clapeyron slope (K/Pa) 
+#define CW    4.22e3   // specific heat capacity of water (J/kg/K)
 
 /*Model processing*/
@@ -176 +178 @@
         /*Get conductivity from inputs*/
         IssmDouble conductivity = GetConductivity(element);
+//if(element->Id()==1){
+//      printf("Conductivity at CreateKMatrix: %g \n",conductivity);
+//}
 
         /* Start  looping on the number of gaussian points: */
@@ -208 +213 @@
         IssmDouble  lr,br,vx,vy,beta;
         IssmDouble  alpha2,frictionheat;
+   IssmDouble  PMPheat,dpressure_water[2],dbed[2];      
         IssmDouble* xyz_list = NULL;
 
@@ -238 +244 @@
         /*Get conductivity from inputs*/
         IssmDouble conductivity = GetConductivity(element);
+//if(element->Id()==1){
+//      printf("Conductivity in CreatePVector: %g \n",conductivity);
+//}
 
         /*Build friction element, needed later: */
@@ -251 +260 @@
                 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);
@@ -283 +293 @@
                 if(pressure_water>pressure_ice) pressure_water = pressure_ice;
 
-                meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
+                /*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=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
+
+                meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(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*
                  (
@@ -315 +330 @@
         IssmDouble* xyz_list = NULL;
 
+        /*Get gravity from parameters*/
+           IssmDouble  g = element->GetMaterialParameter(ConstantsGEnum);
+
         /*Fetch number of nodes for this finite element*/
         int numnodes = element->GetNumberOfNodes();
@@ -323 +341 @@
 
         /*Get thickness and base on nodes to apply cap on water head*/
+   IssmDouble* eff_pressure = xNew<IssmDouble>(numnodes);
         IssmDouble* thickness = xNew<IssmDouble>(numnodes);
         IssmDouble* bed       = xNew<IssmDouble>(numnodes);
@@ -344 +363 @@
                 }
 
+                /*Calculate effective pressure*/
+                eff_pressure[i] = rho_ice*g*thickness[i] - rho_water*g*(values[i]-bed[i]);
+        
                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
                 if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
@@ -350 +372 @@
         /*Add input to the element: */
         element->AddInput(HydrologyHeadEnum,values,element->GetElementType());
+   element->AddInput(EffectivePressureEnum,eff_pressure,P1Enum);
 
         /*Update reynolds number according to new solution*/
@@ -356 +379 @@
         head_input->GetInputDerivativeAverageValue(&dh[0],xyz_list);
         IssmDouble conductivity = GetConductivity(element);
+//if(element->Id()==1){
+//      printf("Conductivity in UpdateInputFromSolution: %g \n",conductivity);
+//}
+
         IssmDouble reynolds = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])/(2.*NU);
         element->AddInput(HydrologyReynoldsEnum,&reynolds,P0Enum);
@@ -365 +392 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<int>(doflist);
+        xDelete<IssmDouble>(eff_pressure);
 }/*}}}*/
 void           HydrologySommersAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
@@ -385 +413 @@
         reynolds_input->GetInputAverage(&reynolds);
         gap_input->GetInputAverage(&gap);
-
+        
         /*Compute conductivity*/
         IssmDouble conductivity = pow(gap,3)*g/(12.*NU*(1+OMEGA*reynolds));
@@ -414 +442 @@
         IssmDouble  alpha2,frictionheat;
         IssmDouble* xyz_list = NULL;
-
+   IssmDouble  dpressure_water[2],dbed[2],PMPheat;
 
         /*Retrieve all inputs and parameters*/
@@ -438 +466 @@
         /*Get conductivity from inputs*/
         IssmDouble conductivity = GetConductivity(element);
-
+//if(element->Id()==1){
+//      printf("Conductivity at gap update: %g \n",conductivity);
+//}
         /*Build friction element, needed later: */
         Friction* friction=new Friction(element,2);
@@ -454 +484 @@
                 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);
@@ -486 +517 @@
                 if(pressure_water>pressure_ice) pressure_water = pressure_ice;
       
-
-                meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
+      /* 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=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
+        
+                meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
                 _assert_(meltrate>0.);
 

issm/trunk-jpl/src/c/cores/hydrology_core.cpp

@@ -88 +88 @@
                 if(save_results){
                         if(VerboseSolution()) _printf0_("   saving results \n");
-                        int outputs[2] = {HydrologyHeadEnum,HydrologyGapHeightEnum};
-                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
+                        int outputs[3] = {HydrologyHeadEnum,HydrologyGapHeightEnum,EffectivePressureEnum};
+                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
                         
                         /*unload results*/