source: issm/oecreview/Archive/16554-17801/ISSM-16599-16600.diff

Last change on this file was 17802, checked in by Mathieu Morlighem, 11 years ago

Added archives
File size: 16.1 KB
RevLine 
[17802]1Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 16599)
4+++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 16600)
5@@ -13,6 +13,7 @@
6 int hydrology_model;
7 int sedimentlimit_flag;
8 int transfer_flag;
9+ int penalty_lock;
10 bool isefficientlayer;
11 IssmDouble sedimentlimit;
12 IssmDouble penalty_factor;
13@@ -30,6 +31,7 @@
14 iomodel->FetchData(&transfer_flag,HydrologydcTransferFlagEnum);
15 iomodel->FetchData(&penalty_factor,HydrologydcPenaltyFactorEnum);
16 iomodel->FetchData(&rel_tol,HydrologydcRelTolEnum);
17+ iomodel->FetchData(&penalty_lock,HydrologydcPenaltyLockEnum);
18
19 if(sedimentlimit_flag==1){
20 iomodel->FetchData(&sedimentlimit,HydrologydcSedimentlimitEnum);
21@@ -47,6 +49,7 @@
22 parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag));
23 parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag));
24 parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol));
25+ parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock));
26
27 }/*}}}*/
28 void HydrologyDCInefficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
29Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
30===================================================================
31--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16599)
32+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16600)
33@@ -120,6 +120,7 @@
34 HydrologydcTransferFlagEnum,
35 HydrologydcLeakageFactorEnum,
36 HydrologydcPenaltyFactorEnum,
37+ HydrologydcPenaltyLockEnum,
38 HydrologyLayerEnum,
39 HydrologySedimentEnum,
40 HydrologyEfficientEnum,
41Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
42===================================================================
43--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16599)
44+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16600)
45@@ -128,6 +128,7 @@
46 case HydrologydcTransferFlagEnum : return "HydrologydcTransferFlag";
47 case HydrologydcLeakageFactorEnum : return "HydrologydcLeakageFactor";
48 case HydrologydcPenaltyFactorEnum : return "HydrologydcPenaltyFactor";
49+ case HydrologydcPenaltyLockEnum : return "HydrologydcPenaltyLock";
50 case HydrologyLayerEnum : return "HydrologyLayer";
51 case HydrologySedimentEnum : return "HydrologySediment";
52 case HydrologyEfficientEnum : return "HydrologyEfficient";
53Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
54===================================================================
55--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16599)
56+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16600)
57@@ -128,6 +128,7 @@
58 else if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum;
59 else if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum;
60 else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum;
61+ else if (strcmp(name,"HydrologydcPenaltyLock")==0) return HydrologydcPenaltyLockEnum;
62 else if (strcmp(name,"HydrologyLayer")==0) return HydrologyLayerEnum;
63 else if (strcmp(name,"HydrologySediment")==0) return HydrologySedimentEnum;
64 else if (strcmp(name,"HydrologyEfficient")==0) return HydrologyEfficientEnum;
65@@ -135,11 +136,11 @@
66 else if (strcmp(name,"WaterTransfer")==0) return WaterTransferEnum;
67 else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum;
68 else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
69- else if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum;
70 else stage=2;
71 }
72 if(stage==2){
73- if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum;
74+ if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum;
75+ else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum;
76 else if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum;
77 else if (strcmp(name,"InversionGradientScaling")==0) return InversionGradientScalingEnum;
78 else if (strcmp(name,"InversionIscontrol")==0) return InversionIscontrolEnum;
79@@ -258,11 +259,11 @@
80 else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
81 else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
82 else if (strcmp(name,"Surface")==0) return SurfaceEnum;
83- else if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum;
84 else stage=3;
85 }
86 if(stage==3){
87- if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum;
88+ if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum;
89+ else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum;
90 else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
91 else if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum;
92 else if (strcmp(name,"SurfaceforcingsS0p")==0) return SurfaceforcingsS0pEnum;
93@@ -381,11 +382,11 @@
94 else if (strcmp(name,"Element")==0) return ElementEnum;
95 else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
96 else if (strcmp(name,"FileParam")==0) return FileParamEnum;
97- else if (strcmp(name,"Input")==0) return InputEnum;
98 else stage=4;
99 }
100 if(stage==4){
101- if (strcmp(name,"IntInput")==0) return IntInputEnum;
102+ if (strcmp(name,"Input")==0) return InputEnum;
103+ else if (strcmp(name,"IntInput")==0) return IntInputEnum;
104 else if (strcmp(name,"InputToExtrude")==0) return InputToExtrudeEnum;
105 else if (strcmp(name,"InputToL2Project")==0) return InputToL2ProjectEnum;
106 else if (strcmp(name,"IntParam")==0) return IntParamEnum;
107@@ -504,11 +505,11 @@
108 else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum;
109 else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
110 else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
111- else if (strcmp(name,"StressTensor")==0) return StressTensorEnum;
112 else stage=5;
113 }
114 if(stage==5){
115- if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
116+ if (strcmp(name,"StressTensor")==0) return StressTensorEnum;
117+ else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
118 else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
119 else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
120 else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
121Index: ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp
122===================================================================
123--- ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 16599)
124+++ ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 16600)
125@@ -639,8 +639,9 @@
126 void Pengrid::ConstraintActivateHydrologyDCInefficient(int* punstable){
127
128 // The penalty is stable if it doesn't change during two consecutive iterations.
129- int unstable = 0;
130+ int unstable=0;
131 int new_active;
132+ int penalty_lock;
133 IssmDouble pressure;
134 IssmDouble h;
135 IssmDouble h_max;
136@@ -655,16 +656,28 @@
137 /*Get sediment water head h*/
138 element->GetInputValue(&h,node,SedimentHeadEnum);
139 element->GetHydrologyDCInefficientHmax(&h_max,node);
140+ parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum);
141+
142 if (h>h_max)
143 new_active=1;
144 else
145 new_active=0;
146
147- if(this->active==new_active)
148- unstable=0;
149- else
150- unstable=1;
151+ if(this->active==new_active){
152+ unstable=0;
153+ }
154+ else{
155+ unstable=1;
156+ if(penalty_lock)zigzag_counter++;
157+ }
158
159+ /*If penalty keeps zigzagging more than penalty_lock times: */
160+ if(penalty_lock){
161+ if(zigzag_counter>penalty_lock){
162+ unstable=0;
163+ active=1;
164+ }
165+ }
166 /*Set penalty flag*/
167 this->active=new_active;
168
169Index: ../trunk-jpl/src/c/classes/FemModel.cpp
170===================================================================
171--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 16599)
172+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 16600)
173@@ -1411,4 +1411,12 @@
174 delete transferg;
175 }
176 /*}}}*/
177+void FemModel::HydrologyEPLThicknessx(void){ /*{{{*/
178+
179+ for (int i=0;i<elements->Size();i++){
180+ Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
181+ element->ComputeEPLThickness();
182+ }
183+}
184+/*}}}*/
185 #endif
186Index: ../trunk-jpl/src/c/classes/Elements/Element.h
187===================================================================
188--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 16599)
189+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 16600)
190@@ -135,6 +135,7 @@
191 virtual void GetHydrologyTransfer(Vector<IssmDouble>* transfer)=0;
192 virtual void HydrologyEPLGetMask(Vector<IssmDouble>* mask)=0;
193 virtual void HydrologyEPLGetActive(Vector<IssmDouble>* active)=0;
194+ virtual void ComputeEPLThickness(void)=0;
195 #endif
196
197 #ifdef _HAVE_GROUNDINGLINE_
198Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
199===================================================================
200--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16599)
201+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16600)
202@@ -7223,6 +7223,12 @@
203 }
204 }
205 /*}}}*/
206+/*FUNCTION Tria::ComputeEPLThickness{{{*/
207+void Tria::ComputeEPLThickness(void){
208+}
209+/*}}}*/
210+
211+
212 #endif
213
214 #ifdef _HAVE_MASSTRANSPORT_
215Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
216===================================================================
217--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16599)
218+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16600)
219@@ -302,6 +302,7 @@
220 void GetHydrologyTransfer(Vector<IssmDouble>* transfer);
221 void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
222 void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
223+ void ComputeEPLThickness(void);
224 bool AllActive(void);
225 bool AnyActive(void);
226 #endif
227Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
228===================================================================
229--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16599)
230+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16600)
231@@ -10844,6 +10844,17 @@
232
233 }
234 /*}}}*/
235+/*FUNCTION Penta::ComputeEPLThickness{{{*/
236+void Penta::ComputeEPLThickness(void){
237+
238+ if (!IsOnBed())return;
239+
240+ Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
241+ tria->ComputeEPLThickness();
242+ delete tria->material; delete tria;
243+
244+}
245+/*}}}*/
246 /*FUNCTION Penta::InputUpdateFromSolutionHydrologyDCInefficient{{{*/
247 void Penta::InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution){
248
249Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
250===================================================================
251--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 16599)
252+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 16600)
253@@ -322,6 +322,7 @@
254 void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
255 void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
256 void InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
257+ void ComputeEPLThickness(void);
258 #endif
259
260 void UpdateConstraintsExtrudeFromBase(void){_error_("not implemented yet");};
261Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
262===================================================================
263--- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 16599)
264+++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 16600)
265@@ -100,6 +100,7 @@
266 void GetHydrologyTransfer(Vector<IssmDouble>* transfer){_error_("not implemented yet");};
267 void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){_error_("not implemented yet");};
268 void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){_error_("not implemented yet");};
269+ void ComputeEPLThickness(void){_error_("not implemented yet");};
270 #endif
271 void GetSolutionFromInputs(Vector<IssmDouble>* solution){_error_("not implemented yet");};
272 void GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum){_error_("not implemented yet");};
273Index: ../trunk-jpl/src/c/classes/FemModel.h
274===================================================================
275--- ../trunk-jpl/src/c/classes/FemModel.h (revision 16599)
276+++ ../trunk-jpl/src/c/classes/FemModel.h (revision 16600)
277@@ -101,6 +101,7 @@
278 #ifdef _HAVE_HYDROLOGY_
279 void HydrologyTransferx(void);
280 void HydrologyEPLupdateDomainx(void);
281+ void HydrologyEPLThicknessx(void);
282 #endif
283 };
284
285Index: ../trunk-jpl/src/m/classes/hydrologydc.m
286===================================================================
287--- ../trunk-jpl/src/m/classes/hydrologydc.m (revision 16599)
288+++ ../trunk-jpl/src/m/classes/hydrologydc.m (revision 16600)
289@@ -8,6 +8,7 @@
290 water_compressibility = 0;
291 isefficientlayer = 0;
292 penalty_factor = 0;
293+ penalty_lock = 0;
294 rel_tol = 0;
295 sedimentlimit_flag = 0;
296 sedimentlimit = 0;
297@@ -106,6 +107,7 @@
298 fielddisplay(obj,'water_compressibility','compressibility of water [Pa^-1]');
299 fielddisplay(obj,'isefficientlayer','do we use an efficient drainage system [1: true; 0: false]');
300 fielddisplay(obj,'penalty_factor','exponent of the value used in the penalisation method [dimensionless]');
301+ fielddisplay(obj,'penalty_lock','stabilize unstable constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
302 fielddisplay(obj,'rel_tol','tolerance of the nonlinear iteration for the transfer between layers [dimensionless]');
303 fielddisplay(obj,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer');
304 disp(sprintf('%55s 0: no limit',' '));
305@@ -146,6 +148,7 @@
306 WriteData(fid,'object',obj,'fieldname','water_compressibility','format','Double');
307 WriteData(fid,'object',obj,'fieldname','isefficientlayer','format','Boolean');
308 WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
309+ WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer');
310 WriteData(fid,'object',obj,'fieldname','rel_tol','format','Double');
311 WriteData(fid,'object',obj,'fieldname','sedimentlimit_flag','format','Integer');
312 WriteData(fid,'object',obj,'fieldname','transfer_flag','format','Integer');
313Index: ../trunk-jpl/src/m/enum/HydrologydcPenaltyLockEnum.m
314===================================================================
315--- ../trunk-jpl/src/m/enum/HydrologydcPenaltyLockEnum.m (revision 0)
316+++ ../trunk-jpl/src/m/enum/HydrologydcPenaltyLockEnum.m (revision 16600)
317@@ -0,0 +1,11 @@
318+function macro=HydrologydcPenaltyLockEnum()
319+%HYDROLOGYDCPENALTYLOCKENUM - Enum of HydrologydcPenaltyLock
320+%
321+% WARNING: DO NOT MODIFY THIS FILE
322+% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
323+% Please read src/c/shared/Enum/README for more information
324+%
325+% Usage:
326+% macro=HydrologydcPenaltyLockEnum()
327+
328+macro=StringToEnum('HydrologydcPenaltyLock');
329Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
330===================================================================
331--- ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16599)
332+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16600)
333@@ -120,6 +120,7 @@
334 def HydrologydcTransferFlagEnum(): return StringToEnum("HydrologydcTransferFlag")[0]
335 def HydrologydcLeakageFactorEnum(): return StringToEnum("HydrologydcLeakageFactor")[0]
336 def HydrologydcPenaltyFactorEnum(): return StringToEnum("HydrologydcPenaltyFactor")[0]
337+def HydrologydcPenaltyLockEnum(): return StringToEnum("HydrologydcPenaltyLock")[0]
338 def HydrologyLayerEnum(): return StringToEnum("HydrologyLayer")[0]
339 def HydrologySedimentEnum(): return StringToEnum("HydrologySediment")[0]
340 def HydrologyEfficientEnum(): return StringToEnum("HydrologyEfficient")[0]
Note: See TracBrowser for help on using the repository browser.