Changeset 19091 for issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

Timestamp:

02/05/15 15:24:10 (10 years ago)

Author:

bdef

Message:

CHG: changes in how we deal with the zigzag counter on the EPL

File:

• issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (modified) (4 diffs)

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

@@ -646 +646 @@
 }/*}}}*/
 
-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;
@@ -672 +672 @@
         int         eplflip_lock;
         int         numnodes      =basalelement->GetNumberOfNodes();
+        IssmDouble  new_active;
+        IssmDouble  recurent;
         IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
         IssmDouble* old_active    =xNew<IssmDouble>(numnodes);
@@ -697 +699 @@
 
         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)*/
@@ -723 +731 @@
                                 /*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;
                 }
         }