Changeset 21208

Timestamp:

09/16/16 09:07:59 (9 years ago)

Author:

aleahsommers

Message:

CHG: Changed creep term on RHS of elliptic eqn to use lagged head from previous time step and added new HydrologyHeadEnum

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologySommersAnalysis.cpp (modified) (8 diffs)
• classes/Elements/Element.cpp (modified) (1 diff)
• cores/hydrology_core.cpp (modified) (1 diff)
• shared/Enum/EnumDefinitions.h (modified) (2 diffs)

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

@@ -120 +120 @@
         iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
         iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
-
+   iomodel->FetchDataToInput(elements,"md.hydrology.head",HydrologyHeadOldEnum); 
+        /*iomodel->FetchDataToInput(elements,"md.hydrology.eff_pressure",EffectivePressureEnum);
+*/ 
         iomodel->FindConstant(&frictionlaw,"md.friction.law");
         /*Friction law variables*/
@@ -210 +212 @@
         /*Intermediaries */
         IssmDouble  Jdet,meltrate,G,dh[2],B,A,n;
-        IssmDouble  gap,bed,thickness,head,ieb;
+        IssmDouble  gap,bed,thickness,head,ieb,head_old;
         IssmDouble  lr,br,vx,vy,beta;
         IssmDouble  alpha2,frictionheat;
@@ -241 +243 @@
         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 conductivity from inputs*/
@@ -270 +273 @@
                 vx_input->GetInputValue(&vx,gauss);
                 vy_input->GetInputValue(&vy,gauss);
+      headold_input->GetInputValue(&head_old,gauss);
 
                 /*Get ice A parameter*/
@@ -293 +297 @@
                 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]);
@@ -304 +312 @@
                  (
                   meltrate*(1/rho_water-1/rho_ice)
-                  +A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap
+                  +A*pow(fabs(pressure_ice - pressure_water_old),n-1)*(pressure_ice - pressure_water_old)*gap
                   -beta*sqrt(vx*vx+vy*vy)
                   +ieb
                   )*basis[i];           
-        }
-
+        }
         /*Clean up and return*/
         xDelete<IssmDouble>(xyz_list);
@@ -379 +386 @@
         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);
@@ -535 +539 @@
         /*Divide by connectivity*/
         newgap = newgap/totalweights;
-        IssmDouble mingap = 0.001;
+        IssmDouble mingap = 1e-5;
         if(newgap<mingap) newgap=mingap;
 

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -1524 +1524 @@
                                 name==HydrologydcEplThicknessEnum ||
                                 name==HydrologydcMaskEplactiveNodeEnum ||
+                                name==HydrologyHeadEnum ||
+                 name==HydrologyHeadOldEnum ||          
                                 name==StressbalanceConvergenceNumStepsEnum || 
                                 name==MeshVertexonbaseEnum 

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

@@ -78 +78 @@
         }
 
-        else if (hydrology_model==HydrologysommersEnum){
-                femmodel->SetCurrentConfiguration(HydrologySommersAnalysisEnum);
-                solutionsequence_nonlinear(femmodel,modify_loads);
+        else if (hydrology_model==HydrologysommersEnum){        
+                femmodel->SetCurrentConfiguration(HydrologySommersAnalysisEnum);        
+      InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum); 
+                solutionsequence_nonlinear(femmodel,modify_loads); 
                 if(VerboseSolution()) _printf0_("   updating gap height\n");
                 HydrologySommersAnalysis* analysis = new HydrologySommersAnalysis();
                 analysis->UpdateGapHeight(femmodel);
-                delete analysis;        
+                delete analysis;
                 
                 if(save_results){

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

@@ -166 +166 @@
         HydrologysommersEnum,
         HydrologyHeadEnum,
+   HydrologyHeadOldEnum,
         HydrologyGapHeightEnum,
         HydrologyBumpSpacingEnum,
@@ -544 +545 @@
         GiaCrossSectionShapeEnum,
         GiadWdtEnum,
-        GiaWEnum,
+        GiaWEnum, 
         /*}}}*/
         /*Results{{{*/