Changeset 23950

Timestamp:

05/29/19 15:35:37 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: changed the way effective pressure is being updated

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyGlaDSAnalysis.cpp (modified) (4 diffs)
• analyses/HydrologyGlaDSAnalysis.h (modified) (1 diff)
• cores/hydrology_core.cpp (modified) (1 diff)

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

@@ -152 +152 @@
 
         /*Get all inputs and parameters*/
-        IssmDouble dt,c_t;
-        element->FindParam(&dt,TimesteppingTimeStepEnum);
-        element->FindParam(&c_t,HydrologyPressureMeltCoefficientEnum);
+        IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
+        IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
+        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
+        IssmDouble g         = element->FindParam(ConstantsGEnum);
         Input* k_input   = element->GetInput(HydrologySheetConductivityEnum);_assert_(k_input);
         Input* phi_input = element->GetInput(HydraulicPotentialEnum);      _assert_(phi_input);
@@ -186 +187 @@
                                                         + coeff*dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
                         }
+                }
+
+                /*Transient term if dt*/
+                if(dt>0.){
                 }
         }
@@ -286 +291 @@
 }/*}}}*/
 void           HydrologyGlaDSAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-
-        /*Intermediary*/
-        int* doflist = NULL;
-        IssmDouble phi_0, phi_m, pi;
-
-        /*Fetch number of nodes for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
-
-        /*Fetch dof list and allocate solution vector*/
-        element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
-        IssmDouble* values = xNew<IssmDouble>(numnodes);
-
-        /*Get thickness and base on nodes to apply cap on water head*/
-   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);
-
-        /*Use the dof list to index into the solution vector: */
-        for(int i=0;i<numnodes;i++){
-                values[i]=solution[doflist[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*/
-                N[i] = phi_0 - values[i];
-
-                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
-                if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
-        }
-
-        /*Add input to the element: */
-        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);
+        element->InputUpdateFromSolutionOneDof(solution,HydraulicPotentialEnum);
+        this->UpdateEffectivePressure(element);
 }/*}}}*/
 void           HydrologyGlaDSAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
@@ -416 +369 @@
         delete gauss;
 }/*}}}*/
+void HydrologyGlaDSAnalysis::UpdateEffectivePressure(FemModel* femmodel){/*{{{*/
+
+        for(int j=0;j<femmodel->elements->Size();j++){
+                Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
+                UpdateEffectivePressure(element);
+        }
+
+}/*}}}*/
+void HydrologyGlaDSAnalysis::UpdateEffectivePressure(Element* element){/*{{{*/
+
+        /*Skip if water or ice shelf element*/
+        if(element->IsFloating()) return;
+        
+        /*Intermediary*/
+        IssmDouble phi_0, phi_m, p_i;
+        IssmDouble H,b,phi;
+
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Get thickness and base on nodes to apply cap on water head*/
+   IssmDouble* N = xNew<IssmDouble>(numnodes);
+        IssmDouble  rho_ice   = element->FindParam(MaterialsRhoIceEnum);
+        IssmDouble  rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
+        IssmDouble  g         = element->FindParam(ConstantsGEnum);
+        Input* H_input   = element->GetInput(ThicknessEnum); _assert_(H_input);
+        Input* b_input   = element->GetInput(BedEnum); _assert_(b_input);
+        Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss();
+        for(int iv=0;iv<numnodes;iv++){
+                gauss->GaussNode(element->FiniteElement(),iv);
+
+                /*Get input values at gauss points*/
+                H_input->GetInputValue(&H,gauss);
+                b_input->GetInputValue(&b,gauss);
+                phi_input->GetInputValue(&phi,gauss);
+
+                /*Elevation potential*/
+                phi_m = rho_water*g*b;
+
+                /*Compute overburden pressure*/
+                p_i = rho_ice*g*H;
+
+                /*Copmute overburden potential*/
+                phi_0 = phi_m + p_i;
+
+                /*Calculate effective pressure*/
+                N[iv] = phi_0 - phi;
+
+                if(xIsNan<IssmDouble>(N[iv])) _error_("NaN found in solution vector");
+                if(xIsInf<IssmDouble>(N[iv])) _error_("Inf found in solution vector");
+        }
+
+        element->AddInput(EffectivePressureEnum,N,element->FiniteElement());
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(N);
+}/*}}}*/

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

@@ -34 +34 @@
                 void UpdateSheetThickness(FemModel* femmodel);
                 void UpdateSheetThickness(Element*  element);
+                void UpdateEffectivePressure(FemModel* femmodel);
+                void UpdateEffectivePressure(Element*  element);
 };
 #endif

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

@@ -193 +193 @@
         /*Using the GlaDS model*/
         else if (hydrology_model==HydrologyGlaDSEnum){
+                HydrologyGlaDSAnalysis* analysis = new HydrologyGlaDSAnalysis();
                 femmodel->SetCurrentConfiguration(HydrologyGlaDSAnalysisEnum);
                 InputDuplicatex(femmodel,HydraulicPotentialEnum,HydraulicPotentialOldEnum);
+                analysis->UpdateEffectivePressure(femmodel);
                 solutionsequence_shakti_nonlinear(femmodel);
                 if(VerboseSolution()) _printf0_("   updating sheet thickness\n");
-                HydrologyGlaDSAnalysis* analysis = new HydrologyGlaDSAnalysis();
                 analysis->UpdateSheetThickness(femmodel);
                 delete analysis;