Changeset 26650

Timestamp:

11/19/21 12:21:28 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: moving EffectivePressure calculation to another function

Location:

issm/trunk-jpl/src/c/analyses

Files:

• HydrologyShaktiAnalysis.cpp (modified) (7 diffs)
• HydrologyShaktiAnalysis.h (modified) (1 diff)

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

@@ -154 +154 @@
    parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.relaxation",HydrologyRelaxationEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.storage",HydrologyStorageEnum));
+
+        /*Deal with friction parameters*/
+        int frictionlaw;
+        iomodel->FindConstant(&frictionlaw,"md.friction.law");
+        if(frictionlaw==6){
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
+        }
+        if(frictionlaw==4){
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
+        }
+        if(frictionlaw==1 || frictionlaw==3 || frictionlaw==7){
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
+        }
+        if(frictionlaw==9){
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
+                parameters->AddObject(new IntParam(FrictionCouplingEnum,0));
+        }
 
   /*Requested outputs*/
@@ -465 +486 @@
 
         /*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);
@@ -496 +516 @@
            values[i] = head_old[i] - relaxation*(head_old[i]-values[i]);
 
-                /*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");
@@ -505 +522 @@
         /*Add input to the element: */
         element->AddInput(HydrologyHeadEnum,values,element->GetElementType());
-   element->AddInput(EffectivePressureEnum,eff_pressure,P1Enum);
 
         /*Update reynolds number according to new solution*/
@@ -519 +535 @@
         IssmDouble reynolds = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])/NU;
         element->AddInput(HydrologyReynoldsEnum,&reynolds,P0Enum);
+
+   /*Compute new effective pressure*/
+   this->UpdateEffectivePressure(element);
 
         /*Free resources:*/
@@ -526 +545 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<int>(doflist);
-        xDelete<IssmDouble>(eff_pressure);
    xDelete<IssmDouble>(head_old);
 }/*}}}*/
@@ -724 +742 @@
         delete gauss;
 }/*}}}*/
+void HydrologyShaktiAnalysis::UpdateEffectivePressure(FemModel* femmodel){/*{{{*/
+
+        for(Object* & object : femmodel->elements->objects){
+                Element* element=xDynamicCast<Element*>(object);
+                UpdateEffectivePressure(element);
+        }
+
+}/*}}}*/
+void HydrologyShaktiAnalysis::UpdateEffectivePressure(Element* element){/*{{{*/
+
+        /*Skip if water or ice shelf element*/
+        if(element->IsAllFloating()) return;
+
+        /*Intermediaries*/
+        IssmDouble bed,thickness,head;
+
+        /* Fetch number of nodes and allocate output*/
+   int numnodes = element->GetNumberOfNodes();
+   IssmDouble* N = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        IssmDouble  g          = element->FindParam(ConstantsGEnum);
+        IssmDouble  rho_ice    = element->FindParam(MaterialsRhoIceEnum);
+        IssmDouble  rho_water  = element->FindParam(MaterialsRhoFreshwaterEnum);
+        Input* head_input      = element->GetInput(HydrologyHeadEnum); _assert_(head_input);
+        Input* thickness_input = element->GetInput(ThicknessEnum);     _assert_(thickness_input);
+        Input* base_input      = element->GetInput(BaseEnum);          _assert_(base_input);
+
+
+   Gauss* gauss=element->NewGauss();
+   for (int i=0;i<numnodes;i++){
+      gauss->GaussNode(element->GetElementType(),i);
+
+                base_input->GetInputValue(&bed,gauss);
+                thickness_input->GetInputValue(&thickness,gauss);
+                head_input->GetInputValue(&head,gauss);
+
+                N[i] = rho_ice*g*thickness - rho_water*g*(head-bed);
+        }
+
+        /*Add new gap as an input*/
+        element->AddInput(EffectivePressureEnum,N,element->GetElementType());
+
+        /*Clean up and return*/
+   xDelete<IssmDouble>(N);
+        delete gauss;
+}/*}}}*/

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

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