Changeset 23661

Timestamp:

01/25/19 03:22:12 (6 years ago)

Author:

bdef

Message:

BUG: bug in the treatment of initialization of inactive nodes

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyDCEfficientAnalysis.cpp (modified) (3 diffs)
• analyses/HydrologyDCInefficientAnalysis.cpp (modified) (5 diffs)
• cores/hydrology_core.cpp (modified) (2 diffs)
• solutionsequences/solutionsequence_hydro_nonlinear.cpp (modified) (8 diffs)

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

@@ -410 +410 @@
 void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
         /*Intermediaries*/
-        bool     active_element;
         int      domaintype;
         Element* basalelement=NULL;
@@ -433 +432 @@
 
         /*Fetch dof list and allocate solution vector*/
-        IssmDouble* eplHeads    = xNew<IssmDouble>(numnodes);
-        IssmDouble* basevalue    = xNew<IssmDouble>(numnodes);
         basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
-
-        Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
-        active_element_input->GetInputValue(&active_element);
-
+        IssmDouble* eplHeads = xNew<IssmDouble>(numnodes);
+        IssmDouble* base     = xNew<IssmDouble>(numnodes);
+
+        basalelement->GetInputListOnVertices(&base[0],BaseEnum);
         /*Use the dof list to index into the solution vector: */
         /*If the EPL is not active we revert to the bedrock elevation*/
-
-        if(active_element){
-                for(int i=0;i<numnodes;i++){
+        for(int i=0;i<numnodes;i++){
+                if(basalelement->nodes[i]->IsActive()){
                         eplHeads[i]=solution[doflist[i]];
-                        if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector");
-                        if(xIsInf<IssmDouble>(eplHeads[i])) _error_("Inf found in solution vector");
-                }
-        }
-        else{
-                //Fixing epl head at bedrock elevation when deactivating
-                basalelement->GetInputListOnVertices(&basevalue[0],BaseEnum);
-                for(int i=0;i<numnodes;i++){
-                        eplHeads[i]=basevalue[i];
-                        if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector");
-                        if(xIsInf<IssmDouble>(eplHeads[i])) _error_("Inf found in solution vector");
-                }
+                }
+                else{
+                        eplHeads[i]=base[i];
+                }
+                if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector");
+                if(xIsInf<IssmDouble>(eplHeads[i])) _error_("Inf found in solution vector");
         }
         /*Add input to the element: */
@@ -463 +453 @@
         /*Free ressources:*/
         xDelete<IssmDouble>(eplHeads);
-        xDelete<IssmDouble>(basevalue);
+        xDelete<IssmDouble>(base);
         xDelete<int>(doflist);
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};

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

@@ -488 +488 @@
 
         /*Intermediaries*/
-        bool     thawed_element;
         int                      domaintype;
         Element* basalelement=NULL;
         bool             converged;
-        int*             doflist        =       NULL;
 
         /*Get basal element*/
@@ -506 +504 @@
                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
         }
+        /*Intermediary*/
+        int* doflist = NULL;
 
         /*Fetch number of nodes for this finite element*/
@@ -515 +515 @@
         IssmDouble* pressure = xNew<IssmDouble>(numnodes);
         IssmDouble* residual = xNew<IssmDouble>(numnodes);
-
-        /*Use the dof list to index into the solution vector: */
-
+        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: */
         for(int i=0;i<numnodes;i++){
-                values[i] =solution[doflist[i]];
+                if(basalelement->nodes[i]->IsActive()){
+                        values[i] =solution[doflist[i]];
+                }
+                else{
+                        values[i] = base[i];
+                }
                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
                 if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
@@ -531 +537 @@
                 IssmDouble  penalty_factor,kmax,kappa,h_max;
                 IssmDouble* thickness = xNew<IssmDouble>(numnodes);
-                IssmDouble* base      = xNew<IssmDouble>(numnodes);
                 IssmDouble* transmitivity = xNew<IssmDouble>(numnodes);
 
@@ -540 +545 @@
                 IssmDouble g              = basalelement->FindParam(ConstantsGEnum);
 
-                basalelement->GetInputListOnVertices(thickness,ThicknessEnum);
-                basalelement->GetInputListOnVertices(base,BaseEnum);
-                basalelement->GetInputListOnVertices(transmitivity,HydrologydcSedimentTransmitivityEnum);
+                basalelement->GetInputListOnVertices(&thickness[0],ThicknessEnum);
+                basalelement->GetInputListOnVertices(&transmitivity[0],HydrologydcSedimentTransmitivityEnum);
+
                 kappa=kmax*pow(10.,penalty_factor);
 
-
-                Input* thawed_element_input = basalelement->GetInput(HydrologydcMaskThawedEltEnum); _assert_(thawed_element_input);
-                thawed_element_input->GetInputValue(&thawed_element);
-
                 for(int i=0;i<numnodes;i++){
-                        /*frozen elements heads are set to base elevation*/
-                        if(!thawed_element){
-                                values[i]=base[i];
-                        }
                         GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->GetNode(i));
                         if(values[i]>h_max) {

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

@@ -64 +64 @@
                 femmodel->HydrologyIDSupdateDomainx(&ThawedNodes);
                 if(ThawedNodes>0){
-                        init_time = time-dt; //getting the time back to the start of the timestep
+                        init_time=time-dt; //getting the time back to the start of the timestep
                         hydrotime=init_time;
                         hydrodt=dt/hydroslices; //computing hydro dt from dt and a divider
@@ -120 +120 @@
                 }
         }
-        else if (hydrology_model==HydrologyshaktiEnum){
+ else if (hydrology_model==HydrologyshaktiEnum){
                 femmodel->SetCurrentConfiguration(HydrologyShaktiAnalysisEnum);
                 InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum);

issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp

@@ -57 +57 @@
         if(!isefficientlayer) hydroconverged=true;
 
-        /*Retrieve inputs as the initial state for the non linear iteration*/
+        /*{{{*//*Retrieve inputs as the initial state for the non linear iteration*/
         GetSolutionFromInputsx(&ug_sed,femmodel);
         Reducevectorgtofx(&uf_sed, ug_sed, femmodel->nodes,femmodel->parameters);
@@ -65 +65 @@
                 femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
                 GetSolutionFromInputsx(&ug_epl,femmodel);
-                /*Initialize the element mask*/
+                /*Initialize the EPL element mask*/
                 inefanalysis->ElementizeEplMask(femmodel);
                 effanalysis->InitZigZagCounter(femmodel);
+                /*Initialize the IDS element mask*/
+                femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
+                inefanalysis->ElementizeIdsMask(femmodel);
         }
+        /*}}}*/
         /*The real computation starts here, outermost loop is on the two layer system*/
         for(;;){
@@ -98 +102 @@
                         /*{{{*//*Loop on the sediment layer to deal with the penalization*/
                         for(;;){
-                                /*{{{*/ /*Core of the computation*/
+                                /*{{{*//*Core of the computation*/
                                 if(VerboseSolution()) _printf0_("Building Sediment Matrix...\n");
                                 SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax,femmodel);
@@ -130 +134 @@
                         sedconverged=false;
 
-                        /*Checking convegence on the value of the sediment head*/
+                        /*Checking convergence on the value of the sediment head*/
                         duf=uf_sed_sub_iter->Duplicate();_assert_(duf);
                         uf_sed_sub_iter->Copy(duf);
@@ -160 +164 @@
                         if(VerboseSolution()) _printf0_("==updating mask...\n");
                         femmodel->HydrologyEPLupdateDomainx(&ThickCount);
-                        inefanalysis->ElementizeEplMask(femmodel);
                         ResetZigzagCounterx(femmodel);
                         InputUpdateFromConstantx(femmodel,false,ConvergedEnum);
@@ -172 +175 @@
                                 femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum);
                                 femmodel->UpdateConstraintsL2ProjectionEPLx(&L2Count);
-                                inefanalysis->ElementizeEplMask(femmodel);
                                 femmodel->parameters->SetParam(EplHeadSlopeXEnum,InputToL2ProjectEnum);
                                 solutionsequence_linear(femmodel);
@@ -182 +184 @@
                                 //updating mask after the computation of the epl thickness (Allow to close too thin EPL)
                                 femmodel->HydrologyEPLupdateDomainx(&ThickCount);
-                                inefanalysis->ElementizeEplMask(femmodel);
                                 /*}}}*/
 
@@ -229 +230 @@
                         }
                 }
-                /*}}}*/ /*End of the global EPL loop*/
-
-                /*{{{*/ /*Now dealing with the convergence of the whole system*/
+                /*}}}*//*End of the global EPL loop*/
+                /*{{{*//*Now dealing with the convergence of the whole system*/
                 if(!hydroconverged){
                         //compute norm(du)/norm(u)