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

Last change on this file since 19102 was 19102, checked in by Mathieu Morlighem, 10 years ago

NEW: added 18296-19100
File size: 7.5 KB
RevLine 
[19102]1Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 19090)
4+++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 19091)
5@@ -645,7 +645,7 @@
6 xDelete<IssmDouble>(dbasis);
7 }/*}}}*/
8
9-void HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, int* eplzigzag_counter, Element* element){
10+void HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, int* eplzigzag_counter, Element* element){
11
12 bool active_element;
13 int i,j;
14@@ -671,6 +671,8 @@
15
16 int eplflip_lock;
17 int numnodes =basalelement->GetNumberOfNodes();
18+ IssmDouble new_active;
19+ IssmDouble recurent;
20 IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
21 IssmDouble* old_active =xNew<IssmDouble>(numnodes);
22 IssmDouble* sedhead =xNew<IssmDouble>(numnodes);
23@@ -696,44 +698,55 @@
24 for(i=1;i<numnodes;i++) if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i];
25
26 for(i=0;i<numnodes;i++){
27- /*Activate EPL if residual is >0 */
28+
29+ /*Dealing with the initial value to define zigzaging, preceding iteration the first time we see the node and current evolving vector after*/
30+ recurence->GetValue(&recurent,basalelement->nodes[i]->Sid());
31+ if (recurent>0)vec_mask->GetValue(&old_active[i],basalelement->nodes[i]->Sid());
32+ new_active=old_active[i];
33+ recurence->SetValue(basalelement->nodes[i]->Sid(),1.,ADD_VAL);
34+
35+ /*Now starting to look at the activations*/
36 if(residual[i]>0.){
37- if(old_active[i]==0) eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
38 vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
39+ new_active=1.0;
40 }
41-
42 /*If mask was already one, keep one or colapse*/
43 else if(old_active[i]>0.){
44-
45 vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
46+ new_active=1.0;
47 /*If epl thickness gets under colapse thickness, close the layer*/
48 if(epl_thickness[i]<colapse_thick){
49- eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
50- /*Avoid flipfloping between open and closed states*/
51- if(eplzigzag_counter[basalelement->nodes[i]->Lid()]<eplflip_lock |eplflip_lock==0){
52- vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
53- epl_thickness[i]=init_thick;
54- }
55+ vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
56+ new_active=0.0;
57 }
58 }
59+
60+ if(basalelement->nodes[i]->Sid()==42){
61+ 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()]);
62+ }
63 /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
64 GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->nodes[i]);
65 if(eplhead[i]>=h_max && active_element){
66 for(j=0;j<numnodes;j++){
67 /*Increase of the domain is on the downstream node in term of sediment head*/
68 if(sedhead[j] == sedheadmin){
69- if(old_active[j]==0) {
70- eplzigzag_counter[basalelement->nodes[j]->Lid()] ++;
71- if(eplzigzag_counter[basalelement->nodes[i]->Lid()]<eplflip_lock |eplflip_lock==0){
72- vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
73- }
74- else{
75- vec_mask->SetValue(basalelement->nodes[j]->Sid(),0.,INS_VAL);
76- }
77- }
78+ vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
79 }
80 }
81 }
82+ /*Now checking for nodes zigzaging between open and close*/
83+ if (old_active[i] != new_active){
84+ eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
85+ /*If node changed too much of state, fix it to it's last known state*/
86+ if(eplzigzag_counter[basalelement->nodes[i]->Lid()]>eplflip_lock & eplflip_lock!=0){
87+ vec_mask->SetValue(basalelement->nodes[i]->Sid(),old_active[i],INS_VAL);
88+ new_active=old_active[i];
89+ }
90+ }
91+ /*If node is now closed bring its thickness back to initial*/
92+ if (new_active==0.){
93+ epl_thickness[i]=init_thick;
94+ }
95 }
96 basalelement->AddInput(HydrologydcEplThicknessEnum,epl_thickness,basalelement->GetElementType());
97
98Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h
99===================================================================
100--- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 19090)
101+++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 19091)
102@@ -39,7 +39,7 @@
103 IssmDouble SedimentStoring(Element* element);
104 IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
105 IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
106- void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, int* eplzigzag_counter, Element* element);
107+ void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, int* eplzigzag_counter, Element* element);
108 void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element);
109 void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode);
110 void ComputeEPLThickness(FemModel* femmodel);
111Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
112===================================================================
113--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 19090)
114+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 19091)
115@@ -91,7 +91,7 @@
116 femmodel->UpdateConstraintsx();
117 femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum);
118
119- /*Reset constraint on the ZigZag Lock, this thing doesn't work, it have to disapear*/
120+ /*Reset constraint on the ZigZag Lock*/
121 ResetConstraintsx(femmodel);
122
123 /*{{{*//*Treating the sediment*/
124Index: ../trunk-jpl/src/c/classes/FemModel.cpp
125===================================================================
126--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 19090)
127+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 19091)
128@@ -1912,6 +1912,7 @@
129 void FemModel::HydrologyEPLupdateDomainx(IssmDouble* pEplcount){ /*{{{*/
130
131 Vector<IssmDouble>* mask = NULL;
132+ Vector<IssmDouble>* recurence = NULL;
133 IssmDouble* serial_mask = NULL;
134 Vector<IssmDouble>* active = NULL;
135 IssmDouble* serial_active = NULL;
136@@ -1921,11 +1922,11 @@
137 HydrologyDCInefficientAnalysis* inefanalysis = new HydrologyDCInefficientAnalysis();
138 /*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/
139 mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
140+ recurence=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
141 this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
142-
143 for (int i=0;i<elements->Size();i++){
144 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
145- effanalysis->HydrologyEPLGetMask(mask,eplzigzag_counter,element);
146+ effanalysis->HydrologyEPLGetMask(mask,recurence,eplzigzag_counter,element);
147 }
148 this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum);
149 /*Assemble and serialize*/
150@@ -1933,6 +1934,7 @@
151 serial_mask=mask->ToMPISerial();
152 xDelete<int>(eplzigzag_counter);
153 delete mask;
154+ delete recurence;
155
156 /*Update Mask*/
157 InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
Note: See TracBrowser for help on using the repository browser.