Changeset 24030

Timestamp:

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

Author:

bdef

Message:

CHG: a better way to reset deactivated hydro nodes

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyDCEfficientAnalysis.cpp (modified) (2 diffs)
• analyses/HydrologyDCInefficientAnalysis.cpp (modified) (8 diffs)
• classes/FemModel.cpp (modified) (5 diffs)
• classes/IoModel.cpp (modified) (1 diff)
• shared/Exp/exp.h (modified) (1 diff)

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

@@ -434 +434 @@
         basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
         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 the EPL is not active we revert to the bedrock elevation when deactivating*/
         for(int i=0;i<numnodes;i++){
-                if(basalelement->nodes[i]->IsActive()){
-                        eplHeads[i]=solution[doflist[i]];
-                }
-                else{
-                        eplHeads[i]=base[i];
-                }
+                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");
@@ -453 +446 @@
         /*Free ressources:*/
         xDelete<IssmDouble>(eplHeads);
-        xDelete<IssmDouble>(base);
         xDelete<int>(doflist);
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};

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){

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -1470 +1470 @@
                 ISSM_MPI_Reduce(&local_icefront_area,&total_icefront_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
                 ISSM_MPI_Bcast(&total_icefront_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-                
+
                 basin_icefront_area[basin-1]=total_icefront_area;
         }
-        
+
         this->parameters->AddObject(new DoubleVecParam(FrontalForcingsBasinIcefrontAreaEnum,basin_icefront_area,numbasins));
         xDelete<IssmDouble>(basin_icefront_area);
-        
+
 }/*}}}*/
 void FemModel::IcefrontMassFluxx(IssmDouble* pM, bool scaled){/*{{{*/
@@ -4808 +4808 @@
         delete active;
 
+
         /*Update node activation accordingly*/
-        int counter =0;
-        for (int i=0;i<nodes->Size();i++){
-                Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
-                if(serial_active[node->Sid()]==1.){
-                        node->Activate();
-                        if(!node->IsClone()) counter++;
-                }
-                else{
-                        node->Deactivate();
-                }
+        int         counter  = 0; //this is probably not acurate but we are only interested in positivity
+        for(int i=0;i<elements->Size();i++){
+                Element    *element  = xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                int         numnodes = element->GetNumberOfNodes();
+                IssmDouble *base     = xNew<IssmDouble>(numnodes);
+                element->GetInputListOnNodes(&base[0],BaseEnum);
+                for(int in=0;in<numnodes;in++){
+                        Node* node=element->GetNode(in);
+                        if(serial_active[node->Sid()]==1.){
+                                node->Activate();
+                                if(!node->IsClone()) counter++;
+                        }
+                        else{
+                                node->Deactivate();
+                                node->ApplyConstraint(0,base[in]);
+                        }
+                }
+                xDelete<IssmDouble>(base);
         }
         xDelete<IssmDouble>(serial_active);
@@ -4828 +4837 @@
         counter=sum_counter;
         *pEplcount = counter;
-        if(VerboseSolution()) _printf0_("   Number of active nodes in EPL layer: "<< counter <<"\n");
+        if(VerboseSolution()) {
+                if(counter==0){
+                        _printf0_("   No nodes are active in EPL layer \n");
+                }
+                else {
+                        _printf0_("   Some active nodes in EPL layer \n");
+                }
+        }
 
         /*Update dof indexings*/
@@ -4881 +4897 @@
 
         /*Update node activation accordingly*/
-        int counter =0;
-        for (int i=0;i<nodes->Size();i++){
-                Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
-                if(serial_active[node->Sid()]==1.){
-                        node->Activate();
-                        if(!node->IsClone()) counter++;
-                }
-                else{
-                        node->Deactivate();
-                }
-        }
-
+        int         counter  = 0; //this is probably not acurate but we are only interested in positivity
+        for(int i=0;i<elements->Size();i++){
+                Element    *element  = xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                int         numnodes = element->GetNumberOfNodes();
+                IssmDouble *base     = xNew<IssmDouble>(numnodes);
+
+                element->GetInputListOnNodes(&base[0],BaseEnum);
+
+                for(int in=0;in<numnodes;in++){
+                        Node* node=element->GetNode(in);
+                        if(serial_active[node->Sid()]==1.){
+                                node->Activate();
+                                if(!node->IsClone()) counter++;
+                        }
+                        else{
+                                node->Deactivate();
+                                node->ApplyConstraint(0,base[in]);
+                        }
+                }
+                xDelete<IssmDouble>(base);
+        }
         xDelete<IssmDouble>(serial_active);
         delete inefanalysis;
@@ -4900 +4925 @@
         counter=sum_counter;
         *pIDScount = counter;
-        if(VerboseSolution()) _printf0_("   Number of active nodes in IDS layer: "<< counter <<"\n");
-
+        if(VerboseSolution()) {
+                if(counter==0){
+                        _printf0_("   No nodes are active in IDS layer \n");
+                }
+                else {
+                        _printf0_("   Some active nodes in IDS layer \n");
+                }
+        }
         /*Update dof indexings*/
         this->UpdateConstraintsx();

issm/trunk-jpl/src/c/classes/IoModel.cpp

@@ -402 +402 @@
                         _printf0_(" Marshalled file is corrupted                                            \n");
                         _printf0_("                                                                         \n");
-                        _printf0_("   * Last record found is : \n");
+                        _printf0_("   * Last record found is :                                              \n");
                         _printf0_("     the corresponding model field has probably been marshalled          \n");
                         _printf0_("     incorrectly                                                         \n");

issm/trunk-jpl/src/c/shared/Exp/exp.h

@@ -80 +80 @@
         /*open domain outline file for reading: */
         if ((fid=fopen(domainname,"r"))==NULL){
-                _error_("could not find file \"" << domainname<<"\". Make sure that the file and path provided exist."); 
+                _error_("could not find file \"" << domainname<<"\". Make sure that the file and path provided exist.");
         }