Changeset 23942

Timestamp:

05/29/19 10:31:40 (6 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added sheet thickness evolution

Location:

issm/trunk-jpl/src/c/analyses

Files:

• HydrologyGlaDSAnalysis.cpp (modified) (1 diff)
• HydrologyGlaDSAnalysis.h (modified) (1 diff)

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

@@ -324 +324 @@
         return;
 }/*}}}*/
+
+/*GlaDS specifics*/
+void HydrologyGlaDSAnalysis::UpdateSheetThickness(FemModel* femmodel){/*{{{*/
+
+        for(int j=0;j<femmodel->elements->Size();j++){
+                Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
+                UpdateSheetThickness(element);
+        }
+
+}/*}}}*/
+void HydrologyGlaDSAnalysis::UpdateSheetThickness(Element* element){/*{{{*/
+
+        /*Skip if water or ice shelf element*/
+        if(element->IsFloating()) return;
+
+        /*Intermediaries */
+        IssmDouble  Jdet,w,v,vx,vy,ub,h_old,N;
+        IssmDouble  A,B,n;
+
+        /*Fetch number vertices for this element*/
+        int numvertices = element->GetNumberOfVertices();
+
+        /*Initialize new sheet thickness*/
+        IssmDouble* h_new = xNew<IssmDouble>(numvertices);
+
+        /*Retrieve all inputs and parameters*/
+        IssmDouble  h_r = element->FindParam(HydrologyBumpHeightEnum);
+        IssmDouble  l_r = element->FindParam(HydrologyBumpSpacingEnum);
+        IssmDouble  dt  = element->FindParam(TimesteppingTimeStepEnum);
+        Input* vx_input = element->GetInput(VxEnum);_assert_(vx_input);
+        Input* vy_input = element->GetInput(VyEnum);_assert_(vy_input);
+        Input*  N_input = element->GetInput(EffectivePressureEnum); _assert_(N_input);
+        Input*  h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
+        Input* B_input  = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
+        Input* n_input  = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss();
+        for(int iv=0;iv<numvertices;iv++){
+                gauss->GaussVertex(iv);
+
+                /*Get input values at gauss points*/
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+                h_input->GetInputValue(&h_old,gauss);
+                B_input->GetInputValue(&B,gauss);
+                n_input->GetInputValue(&n,gauss);
+                N_input->GetInputValue(&N,gauss);
+
+                /*Get basal velocity*/
+                ub = sqrt(vx*vx + vy*vy);
+
+                /*Compute cavity opening w*/
+                w  = 0.;
+                if(h_old<h_r) w = ub*(h_r-h_old)/l_r;
+
+                /*Compute closing rate*/
+                A=pow(B,-n);
+                v = 2./pow(n,n)*A*h_old*pow(fabs(N),n-1.)*N;
+
+                /*Get new sheet thickness*/
+                h_new[iv] += h_old + dt*(w-v);
+        }
+
+        element->AddInput(HydrologySheetThicknessEnum,h_new,P1Enum);
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(h_new);
+        delete gauss;
+}/*}}}*/

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

@@ -30 +30 @@
                 void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
                 void           UpdateConstraints(FemModel* femmodel);
+
+                /*Specific to GlaDS*/
+                void UpdateSheetThickness(FemModel* femmodel);
+                void UpdateSheetThickness(Element*  element);
 };
 #endif