Changeset 24030 for issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

Timestamp:

06/19/19 07:14:34 (6 years ago)

Author:

bdef

Message:

CHG: a better way to reset deactivated hydro nodes

File:

• issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (modified) (8 diffs)

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

@@ -364 +364 @@
         basalelement->FindParam(&smb_model,SmbEnum);
 
-        Input*  sed_head_input                   = basalelement->GetInput(SedimentHeadHydrostepEnum);
-        Input*  epl_head_input                   = basalelement->GetInput(EplHeadHydrostepEnum);
-        Input*  base_input                                       = basalelement->GetInput(BaseEnum);
-        Input*  basal_melt_input                 = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(basal_melt_input);
-        Input*  SedTrans_input                   = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
+        Input* sed_head_input   = basalelement->GetInput(SedimentHeadHydrostepEnum);
+        Input* epl_head_input   = basalelement->GetInput(EplHeadHydrostepEnum);
+        Input* base_input                       = basalelement->GetInput(BaseEnum);
+        Input* basal_melt_input = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(basal_melt_input);
+        Input* SedTrans_input   = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
 
         if(dt!= 0.){
@@ -488 +488 @@
 
         /*Intermediaries*/
-        int                      domaintype;
+        int      domaintype;
         Element* basalelement=NULL;
-        bool             converged;
+        bool     converged;
+        int*     doflist = NULL;
 
         /*Get basal element*/
@@ -504 +505 @@
                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
         }
-        /*Intermediary*/
-        int* doflist = NULL;
-
         /*Fetch number of nodes for this finite element*/
         int numnodes = basalelement->GetNumberOfNodes();
@@ -515 +513 @@
         IssmDouble* pressure = xNew<IssmDouble>(numnodes);
         IssmDouble* residual = xNew<IssmDouble>(numnodes);
-        IssmDouble* base     = xNew<IssmDouble>(numnodes);
-
-        basalelement->GetInputListOnVertices(&base[0],BaseEnum);
-        /*Use the dof list to index into the solution vector, frozen nodes are set to initial value: */
+
+        /*Use the dof list to index into the solution vector reseting to base is done at the deactivate stage: */
         for(int i=0;i<numnodes;i++){
-                if(basalelement->nodes[i]->IsActive()){
-                        values[i] =solution[doflist[i]];
-                }
-                else{
-                        values[i] = base[i];
-                }
+                values[i] =solution[doflist[i]];
                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
                 if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
@@ -537 +528 @@
                 IssmDouble  penalty_factor,kmax,kappa,h_max;
                 IssmDouble* thickness = xNew<IssmDouble>(numnodes);
+                IssmDouble* base = xNew<IssmDouble>(numnodes);
                 IssmDouble* transmitivity = xNew<IssmDouble>(numnodes);
 
@@ -546 +538 @@
 
                 basalelement->GetInputListOnVertices(&thickness[0],ThicknessEnum);
+                basalelement->GetInputListOnVertices(&base[0],BaseEnum);
                 basalelement->GetInputListOnVertices(&transmitivity[0],HydrologydcSedimentTransmitivityEnum);
 
@@ -563 +556 @@
                 }
                 xDelete<IssmDouble>(thickness);
+                xDelete<IssmDouble>(base);
                 xDelete<IssmDouble>(transmitivity);
-                xDelete<IssmDouble>(base);
         }
 
@@ -576 +569 @@
         xDelete<IssmDouble>(residual);
         xDelete<IssmDouble>(pressure);
-        xDelete<IssmDouble>(base);
         xDelete<int>(doflist);
         if(domaintype!=Domain2DhorizontalEnum){