Changeset 27293

Timestamp:

09/28/22 08:43:44 (2 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added hydro speed as output

File:

• issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (modified) (1 diff)

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

@@ -444 +444 @@
 void           HydrologyGlaDSAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
         element->InputUpdateFromSolutionOneDof(solution,HydraulicPotentialEnum);
+
+        /*Compute Hydrology Vx and Vy for time stepping purposes (These inputs do not affect GlaDS)*/
+
+        /*Intermediaries*/
+   IssmDouble  dphi[3],h,k,phi;
+        IssmDouble  oceanLS,iceLS;
+        IssmDouble* xyz_list = NULL;
+
+        /*Hard coded coefficients*/
+        const IssmPDouble alpha = 5./4.;
+        const IssmPDouble beta  = 3./2.;
+
+        /*Fetch number vertices for this element*/
+        int numvertices = element->GetNumberOfVertices();
+
+        /*Initialize new sheet thickness*/
+        IssmDouble* vx = xNew<IssmDouble>(numvertices);
+        IssmDouble* vy = xNew<IssmDouble>(numvertices);
+
+        /*Set to 0 if inactive element*/
+        if(element->IsAllFloating() || !element->IsIceInElement()){
+                for(int iv=0;iv<numvertices;iv++) vx[iv] = 0.;
+                for(int iv=0;iv<numvertices;iv++) vy[iv] = 0.;
+                element->AddInput(HydrologyWaterVxEnum,vx,P1DGEnum);
+                element->AddInput(HydrologyWaterVyEnum,vy,P1DGEnum);
+                xDelete<IssmDouble>(vx);
+                xDelete<IssmDouble>(vy);
+                return;
+        }
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        Input *k_input       = element->GetInput(HydrologySheetConductivityEnum); _assert_(k_input);
+        Input *phi_input     = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
+        Input *h_input       = element->GetInput(HydrologySheetThicknessEnum);    _assert_(h_input);
+        Input *oceanLS_input = element->GetInput(MaskOceanLevelsetEnum);          _assert_(oceanLS_input);
+        Input *iceLS_input   = element->GetInput(MaskIceLevelsetEnum);            _assert_(iceLS_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*/
+      phi_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
+      phi_input->GetInputValue(&phi,gauss);
+      h_input->GetInputValue(&h,gauss);
+      k_input->GetInputValue(&k,gauss);
+                oceanLS_input->GetInputValue(&oceanLS,gauss);
+                iceLS_input->GetInputValue(&iceLS,gauss);
+
+                /*Set sheet thickness to zero if floating or no ice*/
+                if(oceanLS<0. || iceLS>0.){
+                        vx[iv] = 0.;
+         vy[iv] = 0.;
+                }
+                else{
+
+         /*Get norm of gradient of hydraulic potential and make sure it is >0*/
+         IssmDouble normgradphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]);
+         if(normgradphi < AEPS) normgradphi = AEPS;
+
+         IssmDouble coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.)/max(AEPS,h); // divide by h to get speed instead of discharge
+
+                        vx[iv] = -coeff*dphi[0];
+                        vy[iv] = -coeff*dphi[1];
+                }
+        }
+
+        element->AddInput(HydrologyWaterVxEnum,vx,P1DGEnum);
+        element->AddInput(HydrologyWaterVyEnum,vy,P1DGEnum);
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(vx);
+        xDelete<IssmDouble>(vy);
+        delete gauss;
 }/*}}}*/
 void           HydrologyGlaDSAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/