source: issm/oecreview/Archive/25834-26739/ISSM-26462-26463.diff@ 26740

Last change on this file since 26740 was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 15.2 KB
RevLine 
[26740]1Index: ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 26462)
4+++ ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 26463)
5@@ -408,12 +408,9 @@
6 Input* enthalpy_input = element->GetInput(enthalpy_enum); _assert_(enthalpy_input);
7 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
8 Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
9- Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
10- Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
11- Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
12
13 /*Build friction element, needed later: */
14- Friction* friction=new Friction(element,dim);
15+ Friction* friction=new Friction(element,3);
16
17 /******** MELTING RATES ************************************//*{{{*/
18 element->NormalBase(&normal_base[0],xyz_list_base);
19@@ -469,7 +466,7 @@
20
21 /*basal friction*/
22 friction->GetAlpha2(&alpha2,gauss);
23- vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); vz_input->GetInputValue(&vz,gauss);
24+ friction->GetBasalSlidingSpeeds(&vx, &vy, &vz, gauss);
25 basalfriction=alpha2*(vx*vx + vy*vy + vz*vz);
26 geothermalflux_input->GetInputValue(&geothermalflux,gauss);
27 /* -Mb= Fb-(q-q_geo)/((1-w)*L*rho), and (1-w)*rho=rho_ice, cf Aschwanden 2012, eqs.1, 2, 66*/
28@@ -946,9 +943,6 @@
29 element->GetInputValue(&converged,ConvergedEnum);
30 if(dt==0. && !converged) enthalpy_enum=EnthalpyPicardEnum; // use enthalpy from last iteration
31 else enthalpy_enum=EnthalpyEnum; // use enthalpy from last time step
32- Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
33- Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
34- Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
35 Input* enthalpy_input = element->GetInput(enthalpy_enum); _assert_(enthalpy_input);
36 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
37 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
38@@ -987,9 +981,7 @@
39 // although one node might have become temperate. So keep heat flux switched on.
40 geothermalflux_input->GetInputValue(&geothermalflux,gauss);
41 friction->GetAlpha2(&alpha2,gauss);
42- vx_input->GetInputValue(&vx,gauss);
43- vy_input->GetInputValue(&vy,gauss);
44- vz_input->GetInputValue(&vz,gauss);
45+ friction->GetBasalSlidingSpeeds(&vx, &vy, &vz, gauss);
46 basalfriction=alpha2*(vx*vx+vy*vy+vz*vz);
47 heatflux=(basalfriction+geothermalflux)/(rho_ice);
48 scalar=gauss->weight*Jdet*heatflux;
49Index: ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
50===================================================================
51--- ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 26462)
52+++ ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 26463)
53@@ -359,8 +359,6 @@
54 IssmDouble g = element->FindParam(ConstantsGEnum);
55 IssmDouble e_v = element->FindParam(HydrologyEnglacialVoidRatioEnum);
56 Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
57- Input* vx_input = element->GetInput(VxBaseEnum);_assert_(vx_input);
58- Input* vy_input = element->GetInput(VyBaseEnum);_assert_(vy_input);
59 Input* h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
60 Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input);
61 Input* b_input = element->GetInput(BedEnum); _assert_(b_input);
62@@ -382,8 +380,6 @@
63 element->NodalFunctions(basis,gauss);
64
65 /*Get input values at gauss points*/
66- vx_input->GetInputValue(&vx,gauss);
67- vy_input->GetInputValue(&vy,gauss);
68 h_input->GetInputValue(&h,gauss);
69 G_input->GetInputValue(&G,gauss);
70 B_input->GetInputValue(&B,gauss);
71@@ -394,6 +390,7 @@
72 H_input->GetInputValue(&H,gauss);
73
74 /*Get basal velocity*/
75+ friction->GetBasalSlidingSpeeds(&vx, &vy ,gauss);
76 ub = sqrt(vx*vx + vy*vy);
77
78 /*Compute cavity opening w*/
79Index: ../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp
80===================================================================
81--- ../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 26462)
82+++ ../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 26463)
83@@ -118,6 +118,10 @@
84 iomodel->FetchDataToInput(inputs,elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum);
85 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);
86 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum);
87+ if(iomodel->domaintype==Domain2DhorizontalEnum){
88+ iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxBaseEnum);
89+ iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyBaseEnum);
90+ }
91 iomodel->FindConstant(&frictionlaw,"md.friction.law");
92
93 /*Friction law variables*/
94@@ -253,8 +257,6 @@
95 Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
96 Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
97 Input* englacial_input = element->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input);
98- Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
99- Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
100 Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input);
101 Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input);
102 Input* headold_input = element->GetInput(HydrologyHeadOldEnum); _assert_(headold_input);
103@@ -286,8 +288,6 @@
104 englacial_input->GetInputValue(&ieb,gauss);
105 lr_input->GetInputValue(&lr,gauss);
106 br_input->GetInputValue(&br,gauss);
107- vx_input->GetInputValue(&vx,gauss);
108- vy_input->GetInputValue(&vy,gauss);
109 headold_input->GetInputValue(&head_old,gauss);
110
111 /*Get ice A parameter*/
112@@ -303,8 +303,7 @@
113
114 /*Compute frictional heat flux*/
115 friction->GetAlpha2(&alpha2,gauss);
116- vx_input->GetInputValue(&vx,gauss);
117- vy_input->GetInputValue(&vy,gauss);
118+ friction->GetBasalSlidingSpeeds(&vx, &vy, gauss);
119 frictionheat=alpha2*(vx*vx+vy*vy);
120
121 /*Get water and ice pressures*/
122@@ -493,8 +492,6 @@
123 Input* base_input = element->GetInput(BaseEnum); _assert_(base_input);
124 Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
125 Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
126- Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
127- Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
128 Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input);
129 Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input);
130
131@@ -522,8 +519,6 @@
132 head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss);
133 lr_input->GetInputValue(&lr,gauss);
134 br_input->GetInputValue(&br,gauss);
135- vx_input->GetInputValue(&vx,gauss);
136- vy_input->GetInputValue(&vy,gauss);
137
138 /*Get ice A parameter*/
139 B_input->GetInputValue(&B,gauss);
140@@ -538,8 +533,7 @@
141
142 /*Compute frictional heat flux*/
143 friction->GetAlpha2(&alpha2,gauss);
144- vx_input->GetInputValue(&vx,gauss);
145- vy_input->GetInputValue(&vy,gauss);
146+ friction->GetBasalSlidingSpeeds(&vx, &vy, gauss);
147 frictionheat=alpha2*(vx*vx+vy*vy);
148
149 /*Get water and ice pressures*/
150Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
151===================================================================
152--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26462)
153+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26463)
154@@ -3488,7 +3488,7 @@
155
156 /*build friction object, used later on: */
157 /*dim=4 is special for HO, which is actually 2.5D*/
158- Friction* friction=new Friction(element,dim==3?4:1);
159+ Friction* friction=new Friction(element,dim==3?2.5:1);
160
161 /*Recover portion of element that is grounded*/
162 if(!(friction_style==SubelementFriction2Enum)) phi=element->GetGroundedPortion(xyz_list_base);
163Index: ../trunk-jpl/src/c/classes/Loads/Friction.cpp
164===================================================================
165--- ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 26462)
166+++ ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 26463)
167@@ -17,15 +17,75 @@
168 /*Constructors/destructors*/
169 Friction::Friction(){/*{{{*/
170 this->element=NULL;
171- this->dim=0;
172 this->law=0;
173+ this->apply_dim = 1;
174+ this->domaintype=-1;
175+ this->vx_input=NULL;
176+ this->vy_input=NULL;
177+ this->vz_input=NULL;
178 }
179 /*}}}*/
180-Friction::Friction(Element* element_in,int dim_in){/*{{{*/
181+Friction::Friction(Element* element_in){/*{{{*/
182+ /* Determine the dimension according to the domain type automatically.
183+ * There are exceptions, e.g. HO, which needs the user to specify the dimension used in Friciton.*/
184+ this->element=element_in;
185
186+ /* Load necessary parameters */
187+ element_in->FindParam(&this->law,FrictionLawEnum);
188+ element_in->FindParam(&this->domaintype,DomainTypeEnum);
189+
190+ /* Load VxBase and VyBase for this special case */
191+ switch(this->domaintype){
192+ case Domain2DhorizontalEnum:
193+ this->apply_dim = 2;
194+ this->vx_input = element_in->GetInput(VxBaseEnum); _assert_(this->vx_input);
195+ this->vy_input = element_in->GetInput(VyBaseEnum); _assert_(this->vy_input);
196+ this->vz_input = NULL;
197+ break;
198+ case Domain2DverticalEnum:
199+ this->apply_dim = 2;
200+ this->vx_input = element_in->GetInput(VxEnum); _assert_(this->vx_input);
201+ this->vy_input = element_in->GetInput(VyEnum); _assert_(this->vy_input);
202+ this->vz_input = NULL;
203+ break;
204+ case Domain3DEnum:
205+ this->apply_dim = 3;
206+ this->vx_input = element_in->GetInput(VxEnum); _assert_(this->vx_input);
207+ this->vy_input = element_in->GetInput(VyEnum); _assert_(this->vy_input);
208+ this->vz_input = element_in->GetInput(VzEnum); _assert_(this->vz_input);
209+ break;
210+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
211+ }
212+}
213+/*}}}*/
214+Friction::Friction(Element* element_in, double dim){/*{{{*/
215+
216 this->element=element_in;
217- this->dim=dim_in;
218+ this->apply_dim = dim;
219+
220+ /* Load necessary parameters */
221 element_in->FindParam(&this->law,FrictionLawEnum);
222+ element_in->FindParam(&this->domaintype,DomainTypeEnum);
223+
224+ /* Load VxBase and VyBase for this special case */
225+ switch(this->domaintype){
226+ case Domain2DhorizontalEnum:
227+ this->vx_input = element_in->GetInput(VxBaseEnum); _assert_(this->vx_input);
228+ this->vy_input = element_in->GetInput(VyBaseEnum); _assert_(this->vy_input);
229+ this->vz_input = NULL;
230+ break;
231+ case Domain2DverticalEnum:
232+ this->vx_input = element_in->GetInput(VxEnum); _assert_(this->vx_input);
233+ this->vy_input = element_in->GetInput(VyEnum); _assert_(this->vy_input);
234+ this->vz_input = NULL;
235+ break;
236+ case Domain3DEnum:
237+ this->vx_input = element_in->GetInput(VxEnum); _assert_(this->vx_input);
238+ this->vy_input = element_in->GetInput(VyEnum); _assert_(this->vy_input);
239+ this->vz_input = element_in->GetInput(VzEnum); _assert_(this->vz_input);
240+ break;
241+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
242+ }
243 }
244 /*}}}*/
245 Friction::~Friction(){/*{{{*/
246@@ -35,7 +95,7 @@
247 /*methods: */
248 void Friction::Echo(void){/*{{{*/
249 _printf_("Friction:\n");
250- _printf_(" dim: " << this->dim<< "\n");
251+ _printf_(" Domain type: " << this->domaintype<< "\n");
252 }
253 /*}}}*/
254 void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
255@@ -854,33 +914,45 @@
256
257 /*diverse*/
258 IssmDouble vx,vy,vz,vmag;
259+
260+ this->vx_input->GetInputValue(&vx, gauss);
261+ this->vy_input->GetInputValue(&vy, gauss);
262+
263+ if ((this->vz_input == NULL) || (this->apply_dim<3)) vz = 0.0;
264+ else this->vz_input->GetInputValue(&vz, gauss);
265
266- switch(dim){
267- case 1:
268- element->GetInputValue(&vx,gauss,VxEnum);
269- vmag=sqrt(vx*vx);
270- break;
271- case 2:
272- element->GetInputValue(&vx,gauss,VxBaseEnum);
273- element->GetInputValue(&vy,gauss,VyBaseEnum);
274- vmag=sqrt(vx*vx+vy*vy);
275- break;
276- case 3:
277- element->GetInputValue(&vx,gauss,VxEnum);
278- element->GetInputValue(&vy,gauss,VyEnum);
279- element->GetInputValue(&vz,gauss,VzEnum);
280- vmag=sqrt(vx*vx+vy*vy+vz*vz);
281- break;
282- case 4:
283- /* This is for HO 3D case, since it requires the horizontal velocities */
284- element->GetInputValue(&vx,gauss,VxEnum);
285- element->GetInputValue(&vy,gauss,VyEnum);
286- vmag=sqrt(vx*vx+vy*vy);
287- break;
288- default:
289- _error_("not supported");
290- }
291+ if (this->apply_dim<2) vy = 0.0;
292
293+ vmag = sqrt(vx*vx+vy*vy+vz*vz);
294 return vmag;
295+}/*}}}*/
296+void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, Gauss* gauss){/*{{{*/
297
298+ this->vx_input->GetInputValue(pvx, gauss);
299+ /*Checks*/
300+ _assert_(!xIsNan<IssmDouble>(*pvx));
301+ _assert_(!xIsInf<IssmDouble>(*pvx));
302 }/*}}}*/
303+void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, Gauss* gauss){/*{{{*/
304+
305+ this->vx_input->GetInputValue(pvx, gauss);
306+ this->vy_input->GetInputValue(pvy, gauss);
307+ /*Checks*/
308+ _assert_(!xIsNan<IssmDouble>(*pvx));
309+ _assert_(!xIsInf<IssmDouble>(*pvx));
310+ _assert_(!xIsNan<IssmDouble>(*pvy));
311+ _assert_(!xIsInf<IssmDouble>(*pvy));
312+}/*}}}*/
313+void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, IssmDouble* pvz, Gauss* gauss){/*{{{*/
314+
315+ this->vx_input->GetInputValue(pvx, gauss);
316+ this->vy_input->GetInputValue(pvy, gauss);
317+ this->vz_input->GetInputValue(pvz, gauss);
318+ /*Checks*/
319+ _assert_(!xIsNan<IssmDouble>(*pvx));
320+ _assert_(!xIsInf<IssmDouble>(*pvx));
321+ _assert_(!xIsNan<IssmDouble>(*pvy));
322+ _assert_(!xIsInf<IssmDouble>(*pvy));
323+ _assert_(!xIsNan<IssmDouble>(*pvz));
324+ _assert_(!xIsInf<IssmDouble>(*pvz));
325+}/*}}}*/
326Index: ../trunk-jpl/src/c/classes/Loads/Friction.h
327===================================================================
328--- ../trunk-jpl/src/c/classes/Loads/Friction.h (revision 26462)
329+++ ../trunk-jpl/src/c/classes/Loads/Friction.h (revision 26463)
330@@ -14,12 +14,17 @@
331
332 public:
333 Element* element;
334- int dim;
335 int law;
336+ int domaintype;
337+ double apply_dim;
338+ Input* vx_input;
339+ Input* vy_input;
340+ Input* vz_input;
341
342 /*methods: */
343 Friction();
344- Friction(Element* element_in,int dim_in);
345+ Friction(Element* element_in);
346+ Friction(Element* element_in, double dim);
347 ~Friction();
348
349 void Echo(void);
350@@ -46,6 +51,9 @@
351
352 IssmDouble EffectivePressure(Gauss* gauss);
353 IssmDouble VelMag(Gauss* gauss);
354+ void GetBasalSlidingSpeeds(IssmDouble* pvx, Gauss* gauss);
355+ void GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, Gauss* gauss);
356+ void GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, IssmDouble* pvz, Gauss* gauss);
357 };
358
359 #endif /* _FRICTION_H_ */
Note: See TracBrowser for help on using the repository browser.