source: issm/oecreview/Archive/18296-19100/ISSM-19090-19091.diff

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

@@ -645 +645 @@
         xDelete<IssmDouble>(dbasis);
 }/*}}}*/
 
-void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, int* eplzigzag_counter, Element* element){
+void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, int* eplzigzag_counter, Element* element){
 
         bool        active_element;
         int         i,j;
@@ -671 +671 @@
 
         int         eplflip_lock;
         int         numnodes      =basalelement->GetNumberOfNodes();
+        IssmDouble  new_active;
+        IssmDouble  recurent;
         IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
         IssmDouble* old_active    =xNew<IssmDouble>(numnodes);
         IssmDouble* sedhead       =xNew<IssmDouble>(numnodes);
@@ -696 +698 @@
         for(i=1;i<numnodes;i++) if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i];
 
         for(i=0;i<numnodes;i++){
-                /*Activate EPL if residual is >0 */
+
+                /*Dealing with the initial value to define zigzaging, preceding iteration the first time we see the node and current evolving vector after*/
+                recurence->GetValue(&recurent,basalelement->nodes[i]->Sid());
+                if (recurent>0)vec_mask->GetValue(&old_active[i],basalelement->nodes[i]->Sid());
+                new_active=old_active[i];
+                recurence->SetValue(basalelement->nodes[i]->Sid(),1.,ADD_VAL);
+
+                /*Now starting to look at the activations*/
                 if(residual[i]>0.){
-                        if(old_active[i]==0) eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
                         vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
+                        new_active=1.0;
                 }
-
                 /*If mask was already one, keep one or colapse*/
                 else if(old_active[i]>0.){
-                        
                         vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
+                        new_active=1.0;
                         /*If epl thickness gets under colapse thickness, close the layer*/
                         if(epl_thickness[i]<colapse_thick){
-                                eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
-                                /*Avoid flipfloping between open and closed states*/
-                                if(eplzigzag_counter[basalelement->nodes[i]->Lid()]<eplflip_lock |eplflip_lock==0){
-                                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
-                                        epl_thickness[i]=init_thick;
-                                }
+                                vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
+                                new_active=0.0;
                         }
                 }
+
+                if(basalelement->nodes[i]->Sid()==42){
+                        printf("old active is %e, and new active is %e, with counter= %i\n",old_active[i],new_active,eplzigzag_counter[basalelement->nodes[i]->Lid()]);
+                }
                 /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
                 GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->nodes[i]);
                 if(eplhead[i]>=h_max && active_element){
                         for(j=0;j<numnodes;j++){
                                 /*Increase of the domain is on the downstream node in term of sediment head*/
                                 if(sedhead[j] == sedheadmin){
-                                        if(old_active[j]==0) {
-                                                eplzigzag_counter[basalelement->nodes[j]->Lid()] ++;
-                                                if(eplzigzag_counter[basalelement->nodes[i]->Lid()]<eplflip_lock |eplflip_lock==0){
-                                                        vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
-                                                }
-                                                else{
-                                                        vec_mask->SetValue(basalelement->nodes[j]->Sid(),0.,INS_VAL);
-                                                }
-                                        }
+                                        vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
                                 }
                         }
                 }
+                /*Now checking for nodes zigzaging between open and close*/
+                if (old_active[i] != new_active){
+                        eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
+                        /*If node changed too much of state, fix it to it's last known state*/
+                        if(eplzigzag_counter[basalelement->nodes[i]->Lid()]>eplflip_lock & eplflip_lock!=0){
+                                vec_mask->SetValue(basalelement->nodes[i]->Sid(),old_active[i],INS_VAL);
+                                new_active=old_active[i];
+                        }
+                }
+                        /*If node is now closed bring its thickness back to initial*/
+                if (new_active==0.){
+                        epl_thickness[i]=init_thick;
+                }
         }
         basalelement->AddInput(HydrologydcEplThicknessEnum,epl_thickness,basalelement->GetElementType());
 

../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h

@@ -39 +39 @@
                 IssmDouble SedimentStoring(Element* element);
                 IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
                 IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
-                void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, int* eplzigzag_counter, Element* element);
+                void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, int* eplzigzag_counter, Element* element);
                 void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element);
                 void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode);
                 void ComputeEPLThickness(FemModel* femmodel);

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

@@ -91 +91 @@
                 femmodel->UpdateConstraintsx();
                 femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum);
                 
-                /*Reset constraint on the ZigZag Lock, this thing doesn't work, it have to disapear*/
+                /*Reset constraint on the ZigZag Lock*/
                 ResetConstraintsx(femmodel);
 
                 /*{{{*//*Treating the sediment*/

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

@@ -1912 +1912 @@
 void FemModel::HydrologyEPLupdateDomainx(IssmDouble* pEplcount){ /*{{{*/
 
         Vector<IssmDouble>* mask                                                        = NULL;
+        Vector<IssmDouble>* recurence                           = NULL;
         IssmDouble*         serial_mask                         = NULL;
         Vector<IssmDouble>* active                                              = NULL;
         IssmDouble*         serial_active                       = NULL;
@@ -1921 +1922 @@
         HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
         /*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/
         mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
+        recurence=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
         this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); 
-
         for (int i=0;i<elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                effanalysis->HydrologyEPLGetMask(mask,eplzigzag_counter,element);
+                effanalysis->HydrologyEPLGetMask(mask,recurence,eplzigzag_counter,element);
         }
         this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum);
         /*Assemble and serialize*/
@@ -1933 +1934 @@
         serial_mask=mask->ToMPISerial();
         xDelete<int>(eplzigzag_counter);
         delete mask;
+        delete recurence;
 
         /*Update Mask*/
         InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);