source: issm/oecreview/Archive/24684-25833/ISSM-24789-24790.diff

Last change on this file was 25834, checked in by Mathieu Morlighem, 4 years ago

CHG: added 24684-25833
File size: 34.0 KB
RevLine 
[25834]1Index: ../trunk-jpl/src/c/cores/smb_core.cpp
2===================================================================
3--- ../trunk-jpl/src/c/cores/smb_core.cpp (revision 24789)
4+++ ../trunk-jpl/src/c/cores/smb_core.cpp (revision 24790)
5@@ -54,7 +54,7 @@
6
7 /*if yes compute necessary intermediaries and start looping*/
8 if (dtslices>1){
9- int substep;
10+ int substep,smb_averaging;
11 IssmDouble global_time,subtime,yts;
12 IssmDouble dt,subdt;
13
14@@ -61,6 +61,7 @@
15 femmodel->parameters->FindParam(&global_time,TimeEnum);
16 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
17 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
18+ femmodel->parameters->FindParam(&smb_averaging,SmbAveragingEnum);
19
20 subtime=global_time-dt; //getting the time back to the start of the timestep
21 subdt=dt/dtslices; //computing substep from dt and a divider
22@@ -82,7 +83,7 @@
23 }
24 delete analysis;
25 /*averaging the transient input*/
26- femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput);
27+ femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,smb_averaging);
28 /*and reset timesteping variables to original*/
29 femmodel->parameters->SetParam(global_time,TimeEnum);
30 femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
31Index: ../trunk-jpl/src/c/cores/hydrology_core.cpp
32===================================================================
33--- ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 24789)
34+++ ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 24790)
35@@ -62,7 +62,7 @@
36 femmodel->parameters->FindParam(&dtslices,HydrologyStepsPerStepEnum);
37
38 if(dtslices>1){
39- int substep, numaveragedinput;
40+ int substep, numaveragedinput, hydro_averaging;
41 IssmDouble global_time, subtime, yts;
42 IssmDouble dt, subdt;
43
44@@ -69,6 +69,7 @@
45 femmodel->parameters->FindParam(&global_time,TimeEnum);
46 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
47 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
48+ femmodel->parameters->FindParam(&hydro_averaging,HydrologyAveragingEnum);
49
50 subtime=global_time-dt; //getting the time back to the start of the timestep
51 subdt=dt/dtslices; //computing hydro dt from dt and a divider
52@@ -116,7 +117,7 @@
53 femmodel->StackTransientInputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput);
54 }
55 /*averaging the stack*/
56- femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput);
57+ femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,hydro_averaging);
58
59 /*And reseting to global time*/
60 femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
61Index: ../trunk-jpl/src/c/classes/FemModel.h
62===================================================================
63--- ../trunk-jpl/src/c/classes/FemModel.h (revision 24789)
64+++ ../trunk-jpl/src/c/classes/FemModel.h (revision 24790)
65@@ -175,7 +175,7 @@
66 void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count);
67 void InitTransientInputx(int* transientinput_enum,int numoutputs);
68 void StackTransientInputx(int* input_enum,int* transientinput_enum,IssmDouble hydrotime,int numoutputs);
69- void AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs);
70+ void AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs,int averaging_method);
71 void UpdateConstraintsx(void);
72 int UpdateVertexPositionsx(void);
73
74Index: ../trunk-jpl/src/c/classes/Inputs2/TriaInput2.h
75===================================================================
76--- ../trunk-jpl/src/c/classes/Inputs2/TriaInput2.h (revision 24789)
77+++ ../trunk-jpl/src/c/classes/Inputs2/TriaInput2.h (revision 24790)
78@@ -37,7 +37,9 @@
79 TriaInput2* GetTriaInput(){return this;};
80 void GetInputValue(IssmDouble* pvalue,Gauss* gauss);
81 void Scale(IssmDouble scalar);
82+ void Pow(IssmDouble scalar);
83 void AXPY(Input2* xinput,IssmDouble scalar);
84+ void PointWiseMult(Input2* xinput);
85 void Serve(int numindices,int* indices);
86 void Serve(int row,int numindices);
87 void ServeCollapsed(int row,int id0,int in1);
88Index: ../trunk-jpl/src/c/classes/Inputs2/PentaInput2.cpp
89===================================================================
90--- ../trunk-jpl/src/c/classes/Inputs2/PentaInput2.cpp (revision 24789)
91+++ ../trunk-jpl/src/c/classes/Inputs2/PentaInput2.cpp (revision 24790)
92@@ -154,8 +154,8 @@
93 void PentaInput2::SetInput(int interp_in,int row,IssmDouble value_in){/*{{{*/
94
95 _assert_(this);
96- _assert_(row>=0);
97- _assert_(row<this->M);
98+ _assert_(row>=0);
99+ _assert_(row<this->M);
100 _assert_(this->N==1);
101
102 this->values[row] = value_in;
103@@ -169,8 +169,8 @@
104 _assert_(this->N==1);
105 for(int i=0;i<numindices;i++){
106 int row = indices[i];
107- _assert_(row>=0);
108- _assert_(row<this->M);
109+ _assert_(row>=0);
110+ _assert_(row<this->M);
111 this->values[row] = values_in[i];
112 }
113 }
114@@ -178,8 +178,8 @@
115 this->Reset(interp_in);
116 for(int i=0;i<numindices;i++){
117 int row = indices[i];
118- _assert_(row>=0);
119- _assert_(row<this->M);
120+ _assert_(row>=0);
121+ _assert_(row<this->M);
122 this->values[row] = values_in[i];
123 }
124 }
125@@ -212,8 +212,8 @@
126
127 for(int i=0;i<numindices;i++){
128 int row = indices[i];
129- _assert_(row>=0);
130- _assert_(row<this->M);
131+ _assert_(row>=0);
132+ _assert_(row<this->M);
133 this->element_values[i] = this->values[row];
134 }
135
136@@ -388,6 +388,12 @@
137 for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*this->element_values[i];
138 }
139 /*}}}*/
140+void PentaInput2::Pow(IssmDouble alpha){/*{{{*/
141+
142+ for(int i=0;i<this->M*this->N;i++) this->values[i] = pow(this->values[i],alpha);
143+ for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = pow(this->element_values[i],alpha);
144+}
145+/*}}}*/
146 void PentaInput2::AXPY(Input2* xinput,IssmDouble alpha){/*{{{*/
147
148 /*xinput is of the same type, so cast it: */
149@@ -400,5 +406,20 @@
150 for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*xpentainput->element_values[i] + this->element_values[i];
151 }
152 /*}}}*/
153+void PentaInput2::PointWiseMult(Input2* xinput){/*{{{*/
154
155+ /*xinput is of the same type, so cast it: */
156+ if(xinput->ObjectEnum()!=PentaInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
157+ PentaInput2* xpentainput=xDynamicCast<PentaInput2*>(xinput);
158+ if(xpentainput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
159+
160+ /* we need to check that the vector sizes are identical*/
161+ if(xpentainput->M!=this->M||xpentainput->N!=this->N) _error_("Operation not permitted because the inputs have different sizes");
162+
163+ /*Carry out the pointwise operation depending on type:*/
164+ for(int i=0;i<this->M*this->N;i++) this->values[i] = xpentainput->values[i] * this->values[i];
165+ for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = xpentainput->element_values[i] * this->element_values[i];
166+}
167+/*}}}*/
168+
169 /*Object functions*/
170Index: ../trunk-jpl/src/c/classes/Inputs2/Inputs2.h
171===================================================================
172--- ../trunk-jpl/src/c/classes/Inputs2/Inputs2.h (revision 24789)
173+++ ../trunk-jpl/src/c/classes/Inputs2/Inputs2.h (revision 24790)
174@@ -53,7 +53,7 @@
175 SegInput2* GetSegInput(int enum_type);
176 TriaInput2* GetTriaInput(int enum_type);
177 TriaInput2* GetTriaInput(int enum_type,IssmDouble time);
178- TriaInput2* GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time);
179+ TriaInput2* GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time,int averaging_method);
180 PentaInput2* GetPentaInput(int enum_type);
181 PentaInput2* GetPentaInput(int enum_type,IssmDouble time);
182 TransientInput2* GetTransientInput(int enum_type);
183Index: ../trunk-jpl/src/c/classes/Inputs2/PentaInput2.h
184===================================================================
185--- ../trunk-jpl/src/c/classes/Inputs2/PentaInput2.h (revision 24789)
186+++ ../trunk-jpl/src/c/classes/Inputs2/PentaInput2.h (revision 24790)
187@@ -36,7 +36,9 @@
188 PentaInput2* GetPentaInput(){return this;};
189 void GetInputValue(IssmDouble* pvalue,Gauss* gauss);
190 void Scale(IssmDouble scalar);
191+ void Pow(IssmDouble scalar);
192 void AXPY(Input2* xinput,IssmDouble scalar);
193+ void PointWiseMult(Input2* xinput);
194 void Serve(int numindices,int* indices);
195 void Serve(int row,int numindices);
196 void ServeCollapsed(int row,int state);
197Index: ../trunk-jpl/src/c/classes/Inputs2/TransientInput2.h
198===================================================================
199--- ../trunk-jpl/src/c/classes/Inputs2/TransientInput2.h (revision 24789)
200+++ ../trunk-jpl/src/c/classes/Inputs2/TransientInput2.h (revision 24790)
201@@ -49,7 +49,7 @@
202 void GetAllTimes(IssmDouble** ptimesteps,int* pnumtimesteps);
203 TriaInput2* GetTriaInput();
204 TriaInput2* GetTriaInput(IssmDouble time);
205- TriaInput2* GetTriaInput(IssmDouble start_time,IssmDouble end_time);
206+ TriaInput2* GetTriaInput(IssmDouble start_time,IssmDouble end_time,int averaging_method);
207 TriaInput2* GetTriaInput(int offset);
208 PentaInput2* GetPentaInput();
209 PentaInput2* GetPentaInput(IssmDouble time);
210@@ -58,7 +58,7 @@
211 IssmDouble GetTimeByOffset(int offset);
212 int GetTimeInputOffset(IssmDouble time);
213 void SetCurrentTimeInput(IssmDouble time);
214- void SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time);
215+ void SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time,int averaging_method);
216 /*numerics:*/
217
218 };
219Index: ../trunk-jpl/src/c/classes/Inputs2/SegInput2.cpp
220===================================================================
221--- ../trunk-jpl/src/c/classes/Inputs2/SegInput2.cpp (revision 24789)
222+++ ../trunk-jpl/src/c/classes/Inputs2/SegInput2.cpp (revision 24790)
223@@ -136,8 +136,8 @@
224 void SegInput2::SetInput(int interp_in,int row,IssmDouble value_in){/*{{{*/
225
226 _assert_(this);
227- _assert_(row>=0);
228- _assert_(row<this->M);
229+ _assert_(row>=0);
230+ _assert_(row<this->M);
231 _assert_(this->N==1);
232
233 this->values[row] = value_in;
234@@ -151,8 +151,8 @@
235 _assert_(this->N==1);
236 for(int i=0;i<numindices;i++){
237 int row = indices[i];
238- _assert_(row>=0);
239- _assert_(row<this->M);
240+ _assert_(row>=0);
241+ _assert_(row<this->M);
242 this->values[row] = values_in[i];
243 }
244 }
245@@ -160,8 +160,8 @@
246 _assert_(this->N==1);
247 for(int i=0;i<numindices;i++){
248 int row = indices[i];
249- _assert_(row>=0);
250- _assert_(row<this->M);
251+ _assert_(row>=0);
252+ _assert_(row<this->M);
253 this->values[row] = values_in[i];
254 }
255 }
256@@ -169,8 +169,8 @@
257 this->Reset(interp_in);
258 for(int i=0;i<numindices;i++){
259 int row = indices[i];
260- _assert_(row>=0);
261- _assert_(row<this->M);
262+ _assert_(row>=0);
263+ _assert_(row<this->M);
264 this->values[row] = values_in[i];
265 }
266 }
267@@ -201,8 +201,8 @@
268
269 for(int i=0;i<numindices;i++){
270 int row = indices[i];
271- _assert_(row>=0);
272- _assert_(row<this->M);
273+ _assert_(row>=0);
274+ _assert_(row<this->M);
275 this->element_values[i] = this->values[row];
276 }
277
278@@ -308,17 +308,38 @@
279 for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*this->element_values[i];
280 }
281 /*}}}*/
282+void SegInput2::Pow(IssmDouble alpha){/*{{{*/
283+
284+ for(int i=0;i<this->M*this->N;i++) this->values[i] = pow(this->values[i],alpha);
285+ for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = pow(this->element_values[i],alpha);
286+}
287+/*}}}*/
288 void SegInput2::AXPY(Input2* xinput,IssmDouble alpha){/*{{{*/
289
290 /*xinput is of the same type, so cast it: */
291 if(xinput->ObjectEnum()!=SegInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
292- SegInput2* xtriainput=xDynamicCast<SegInput2*>(xinput);
293- if(xtriainput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
294+ SegInput2* xseginput=xDynamicCast<SegInput2*>(xinput);
295+ if(xseginput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
296
297 /*Carry out the AXPY operation depending on type:*/
298- for(int i=0;i<this->M*this->N;i++) this->values[i] = alpha*xtriainput->values[i] + this->values[i];
299- for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*xtriainput->element_values[i] + this->element_values[i];
300+ for(int i=0;i<this->M*this->N;i++) this->values[i] = alpha*xseginput->values[i] + this->values[i];
301+ for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*xseginput->element_values[i] + this->element_values[i];
302 }
303 /*}}}*/
304+void SegInput2::PointWiseMult(Input2* xinput){/*{{{*/
305
306+ /*xinput is of the same type, so cast it: */
307+ if(xinput->ObjectEnum()!=SegInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
308+ SegInput2* xseginput=xDynamicCast<SegInput2*>(xinput);
309+ if(xseginput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
310+
311+ /* we need to check that the vector sizes are identical*/
312+ if(xseginput->M!=this->M||xseginput->N!=this->N) _error_("Operation not permitted because the inputs have different sizes");
313+
314+ /*Carry out the AXPY operation depending on type:*/
315+ for(int i=0;i<this->M*this->N;i++) this->values[i] = xseginput->values[i] * this->values[i];
316+ for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = xseginput->element_values[i] * this->element_values[i];
317+}
318+/*}}}*/
319+
320 /*Object functions*/
321Index: ../trunk-jpl/src/c/classes/Inputs2/SegInput2.h
322===================================================================
323--- ../trunk-jpl/src/c/classes/Inputs2/SegInput2.h (revision 24789)
324+++ ../trunk-jpl/src/c/classes/Inputs2/SegInput2.h (revision 24790)
325@@ -34,7 +34,9 @@
326 SegInput2* GetSegInput(){return this;};
327 void GetInputValue(IssmDouble* pvalue,Gauss* gauss);
328 void Scale(IssmDouble scalar);
329+ void Pow(IssmDouble scalar);
330 void AXPY(Input2* xinput,IssmDouble scalar);
331+ void PointWiseMult(Input2* xinput);
332 void Serve(int numindices,int* indices);
333 void Serve(int row,int numindices);
334 int GetResultArraySize(void);
335Index: ../trunk-jpl/src/c/classes/Inputs2/TriaInput2.cpp
336===================================================================
337--- ../trunk-jpl/src/c/classes/Inputs2/TriaInput2.cpp (revision 24789)
338+++ ../trunk-jpl/src/c/classes/Inputs2/TriaInput2.cpp (revision 24790)
339@@ -141,8 +141,8 @@
340 void TriaInput2::SetInput(int interp_in,int row,IssmDouble value_in){/*{{{*/
341
342 _assert_(this);
343- _assert_(row>=0);
344- _assert_(row<this->M);
345+ _assert_(row>=0);
346+ _assert_(row<this->M);
347 _assert_(this->N==1);
348
349 this->values[row] = value_in;
350@@ -156,8 +156,8 @@
351 _assert_(this->N==1);
352 for(int i=0;i<numindices;i++){
353 int row = indices[i];
354- _assert_(row>=0);
355- _assert_(row<this->M);
356+ _assert_(row>=0);
357+ _assert_(row<this->M);
358 this->values[row] = values_in[i];
359 }
360 }
361@@ -165,8 +165,8 @@
362 _assert_(this->N==1);
363 for(int i=0;i<numindices;i++){
364 int row = indices[i];
365- _assert_(row>=0);
366- _assert_(row<this->M);
367+ _assert_(row>=0);
368+ _assert_(row<this->M);
369 this->values[row] = values_in[i];
370 }
371 }
372@@ -174,8 +174,8 @@
373 this->Reset(interp_in);
374 for(int i=0;i<numindices;i++){
375 int row = indices[i];
376- _assert_(row>=0);
377- _assert_(row<this->M);
378+ _assert_(row>=0);
379+ _assert_(row<this->M);
380 this->values[row] = values_in[i];
381 }
382 }
383@@ -206,8 +206,8 @@
384
385 for(int i=0;i<numindices;i++){
386 int row = indices[i];
387- _assert_(row>=0);
388- _assert_(row<this->M);
389+ _assert_(row>=0);
390+ _assert_(row<this->M);
391 this->element_values[i] = this->values[row];
392 }
393
394@@ -363,6 +363,12 @@
395 for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*this->element_values[i];
396 }
397 /*}}}*/
398+void TriaInput2::Pow(IssmDouble alpha){/*{{{*/
399+
400+ for(int i=0;i<this->M*this->N;i++) this->values[i] = pow(this->values[i],alpha);
401+ for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = pow(this->element_values[i],alpha);
402+}
403+/*}}}*/
404 void TriaInput2::AXPY(Input2* xinput,IssmDouble alpha){/*{{{*/
405
406 /*xinput is of the same type, so cast it: */
407@@ -375,5 +381,20 @@
408 for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*xtriainput->element_values[i] + this->element_values[i];
409 }
410 /*}}}*/
411+void TriaInput2::PointWiseMult(Input2* xinput){/*{{{*/
412
413+ /*xinput is of the same type, so cast it: */
414+ if(xinput->ObjectEnum()!=TriaInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
415+ TriaInput2* xtriainput=xDynamicCast<TriaInput2*>(xinput);
416+ if(xtriainput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
417+
418+ /* we need to check that the vector sizes are identical*/
419+ if(xtriainput->M!=this->M||xtriainput->N!=this->N) _error_("Operation not permitted because the inputs have different sizes");
420+
421+ /*Carry out the AXPY operation depending on type:*/
422+ for(int i=0;i<this->M*this->N;i++) this->values[i] = xtriainput->values[i] * this->values[i];
423+ for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = xtriainput->element_values[i] * this->element_values[i];
424+}
425+/*}}}*/
426+
427 /*Object functions*/
428Index: ../trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp
429===================================================================
430--- ../trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp (revision 24789)
431+++ ../trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp (revision 24790)
432@@ -343,7 +343,7 @@
433 return input->GetTriaInput();
434 }
435 }/*}}}*/
436-TriaInput2* Inputs2::GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time){/*{{{*/
437+TriaInput2* Inputs2::GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time,int averaging_method){/*{{{*/
438
439 /*Get input id*/
440 int id = EnumToIndex(enum_in);
441@@ -353,7 +353,7 @@
442 if(!input) return NULL;
443
444 if(input->ObjectEnum()==TransientInput2Enum){
445- return xDynamicCast<TransientInput2*>(input)->GetTriaInput(start_time,end_time);
446+ return xDynamicCast<TransientInput2*>(input)->GetTriaInput(start_time,end_time,averaging_method);
447 }
448 else{
449 _error_("Input "<<EnumToStringx(enum_in)<<" is not an TransientInput2");
450Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
451===================================================================
452--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 24789)
453+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 24790)
454@@ -177,7 +177,7 @@
455 void AddInput2(int input_enum, IssmDouble* values, int interpolation_enum);
456 void AddControlInput(int input_enum,Inputs2* inputs2,IoModel* iomodel,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id);
457 void DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,Inputs2* inputs2,IoModel* iomodel,int input_enum);
458- void CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble init_time,IssmDouble end_time);
459+ void CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int averaging_method);
460 void GetInputAveragesUpToCurrentTime(int input_enum,IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime);
461 IssmDouble GetArea(void);
462 IssmDouble GetHorizontalSurfaceArea(void);
463@@ -188,7 +188,7 @@
464 int GetElementType(void);
465 Input2* GetInput2(int enumtype);
466 Input2* GetInput2(int enumtype,IssmDouble time);
467- Input2* GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time);
468+ Input2* GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method);
469 DatasetInput2* GetDatasetInput2(int inputenum);
470 void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
471 void GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype);
472Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
473===================================================================
474--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 24789)
475+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 24790)
476@@ -1007,7 +1007,7 @@
477 this->AddInput2(distanceenum,&ls[0],P1Enum);
478 }
479 /*}}}*/
480-void Penta::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time){/*{{{*/
481+void Penta::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time,int averaging_method){/*{{{*/
482
483 _assert_(end_time>start_time);
484
485@@ -1016,9 +1016,7 @@
486 IssmDouble current_values[NUMVERTICES];
487 IssmDouble dt;
488 int found,start_offset,end_offset;
489- int averaging_method=0;
490
491-
492 /*Get transient input time steps*/
493 int numtimesteps;
494 IssmDouble *timesteps = NULL;
495Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
496===================================================================
497--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 24789)
498+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 24790)
499@@ -67,7 +67,7 @@
500 void ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
501 void CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum);
502 ElementMatrix* CreateBasalMassMatrix(void);
503- void CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time);
504+ void CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time,int averaging_method);
505 void ElementResponse(IssmDouble* presponse,int response_enum);
506 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
507 int FiniteElement(void);
508@@ -85,7 +85,7 @@
509 void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
510 Input2* GetInput2(int enumtype);
511 Input2* GetInput2(int enumtype,IssmDouble time);
512- Input2* GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");};
513+ Input2* GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method){_error_("not implemented yet!");};
514 void GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
515 void GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
516 DatasetInput2* GetDatasetInput2(int inputenum);
517Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
518===================================================================
519--- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 24789)
520+++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 24790)
521@@ -64,7 +64,7 @@
522 void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
523 Input2* GetInput2(int enumtype);
524 Input2* GetInput2(int enumtype,IssmDouble time);
525- Input2* GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");};
526+ Input2* GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method){_error_("not implemented yet!");};
527 void GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
528 void GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
529 void GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){_error_("not implemented yet");};
530Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h
531===================================================================
532--- ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24789)
533+++ ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24790)
534@@ -68,7 +68,7 @@
535 void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
536 Input2* GetInput2(int enumtype);
537 Input2* GetInput2(int enumtype,IssmDouble time);
538- Input2* GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");};
539+ Input2* GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method){_error_("not implemented yet!");};
540 void GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
541 void GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
542 void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
543Index: ../trunk-jpl/src/c/classes/Elements/Element.h
544===================================================================
545--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 24789)
546+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 24790)
547@@ -233,7 +233,7 @@
548 virtual void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0;
549 virtual void ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum)=0;
550 virtual void CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum){_error_("not implemented yet");};
551- virtual void CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble init_time,IssmDouble end_time){_error_("not implemented yet "<<this->ObjectEnum());};
552+ virtual void CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int averaging_method){_error_("not implemented yet "<<this->ObjectEnum());};
553 virtual void ElementResponse(IssmDouble* presponse,int response_enum)=0;
554 virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
555 virtual int FiniteElement(void)=0;
556@@ -249,7 +249,7 @@
557 virtual DatasetInput2* GetDatasetInput2(int inputenum){_error_("not implemented");};
558 virtual Input2* GetInput2(int inputenum)=0;
559 virtual Input2* GetInput2(int inputenum,IssmDouble time)=0;
560- virtual Input2* GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time)=0;
561+ virtual Input2* GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method)=0;
562 virtual void GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){_error_("not implemented yet");};
563 virtual void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
564 virtual void GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value)=0;
565Index: ../trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp
566===================================================================
567--- ../trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp (revision 24789)
568+++ ../trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp (revision 24790)
569@@ -289,10 +289,10 @@
570
571 }
572 /*}}}*/
573-TriaInput2* TransientInput2::GetTriaInput(IssmDouble start_time, IssmDouble end_time){/*{{{*/
574+TriaInput2* TransientInput2::GetTriaInput(IssmDouble start_time, IssmDouble end_time, int averaging_method){/*{{{*/
575
576 /*Set current time input*/
577- this->SetAverageAsCurrentTimeInput(start_time,end_time);
578+ this->SetAverageAsCurrentTimeInput(start_time,end_time,averaging_method);
579 _assert_(this->current_input);
580
581 /*Cast and return*/
582@@ -420,13 +420,12 @@
583 }
584
585 }/*}}}*/
586-void TransientInput2::SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time){/*{{{*/
587+void TransientInput2::SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time, int averaging_method){/*{{{*/
588
589- IssmDouble dt;
590+ IssmDouble dt,durinv;
591 IssmPDouble eps=1.0e-6;
592 IssmDouble dtsum=0;
593 int found,start_offset,end_offset;
594- int averaging_method = 0;
595
596 /*go through the timesteps, and grab offset for start and end*/
597 IssmDouble temp = start_time-eps;
598@@ -462,9 +461,25 @@
599 }
600 break;
601 case 1: /*Geometric mean*/
602- _error_("Geometric not implemented yet");
603+ if(offset==start_offset){
604+ this->current_input = stepinput->copy();
605+ this->current_input->Scale(dt);
606+ }
607+ else{
608+ this->stepinput->Scale(dt);
609+ this->current_input->PointWiseMult(stepinput)
610+ }
611+ break;
612 case 2: /*Harmonic mean*/
613- _error_("Harmonic not implemented yet");
614+ if(offset==start_offset){
615+ this->current_input = stepinput->copy();
616+ this->current_input->Pow(-1);
617+ this->current_input->Scale(dt);
618+ }
619+ else{
620+ this->stepinput->Pow(-1);
621+ this->current_input->AXPY(stepinput,dt);
622+ }
623 default:
624 _error_("averaging method is not recognised");
625 }
626@@ -471,15 +486,19 @@
627 dtsum+=dt;
628 offset+=1;
629 }
630- /*Integration done, now normalize*/
631+ _assert_(dtsum>0);
632+ durinv=1./dtsum;
633+ /*Integration done, now normalize*/
634 switch(averaging_method){
635 case 0: //Arithmetic mean
636- this->current_input->Scale(1/(dtsum));
637+ this->current_input->Scale(durinv);
638 break;
639 case 1: /*Geometric mean*/
640- _error_("Geometric not implemented yet");
641+ this->current_input->Pow(durinv);
642+ break;
643 case 2: /*Harmonic mean*/
644- _error_("Harmonic not implemented yet");
645+ this->current_input->Scale(durinv);
646+ this->current_input->Pow(-1);
647 default:
648 _error_("averaging method is not recognised");
649 }
650Index: ../trunk-jpl/src/c/classes/FemModel.cpp
651===================================================================
652--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 24789)
653+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 24790)
654@@ -5274,12 +5274,12 @@
655 }
656 }
657 /*}}}*/
658-void FemModel::AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs){ /*{{{*/
659+void FemModel::AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs,int averaging_method){ /*{{{*/
660
661 for(int i=0;i<numoutputs;i++){
662 for(int j=0;j<this->elements->Size();j++){
663 Element* element = xDynamicCast<Element*>(elements->GetObjectByOffset(j));
664- element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time);
665+ element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time,averaging_method);
666 }
667 }
668 }/*}}}*/
669Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
670===================================================================
671--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24789)
672+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24790)
673@@ -1912,7 +1912,7 @@
674 return input;
675 }
676 }/*}}}*/
677-Input2* Tria::GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time){/*{{{*/
678+Input2* Tria::GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time, int averaging_method){/*{{{*/
679
680 /*Get Input from dataset*/
681 if(this->iscollapsed){
682@@ -1919,7 +1919,7 @@
683 _error_("Get Average input not implemented in Penta yet");
684 }
685 else{
686- TriaInput2* input = this->inputs2->GetTriaInput(inputenum,start_time, end_time);
687+ TriaInput2* input = this->inputs2->GetTriaInput(inputenum,start_time, end_time, int averaging_method);
688 if(!input) return input;
689
690 this->InputServe(input);
691@@ -2108,7 +2108,7 @@
692
693 return datasetinput;
694 }/*}}}*/
695-void Tria::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time){/*{{{*/
696+void Tria::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time, int averaging_method){/*{{{*/
697
698 _assert_(end_time>start_time);
699
700@@ -2117,9 +2117,7 @@
701 IssmDouble current_values[NUMVERTICES];
702 IssmDouble dt;
703 int found,start_offset,end_offset;
704- int averaging_method=0;
705
706-
707 /*Get transient input time steps*/
708 int numtimesteps;
709 IssmDouble *timesteps = NULL;
710@@ -2212,7 +2210,7 @@
711 for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = pow(averaged_values[iv], 1./(end_time - start_time));
712 break;
713 case 2: /*Harmonic mean*/
714- for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = 1./(averaged_values[iv]/(end_time - start_time));
715+ for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = (end_time - start_time)/(averaged_values[iv];
716 break;
717 default:
718 _error_("averaging method is not recognised");
719Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
720===================================================================
721--- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 24789)
722+++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 24790)
723@@ -335,6 +335,7 @@
724 int smb_model;
725 int smbsubstepping;
726 int hydrologysubstepping;
727+ int smb_averaging;
728 IssmDouble dt,scalar,sediment_storing;
729 IssmDouble water_head,sediment_transmitivity;
730 IssmDouble water_load,runoff_value,transfer;
731@@ -358,6 +359,7 @@
732 basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
733 basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
734 basalelement->FindParam(&smb_model,SmbEnum);
735+ basalelement->FindParam(&smb_averaging,SmbAveragingEnum);
736
737 Input2* sed_head_input = basalelement->GetInput2(SedimentHeadSubstepEnum);
738 Input2* epl_head_input = basalelement->GetInput2(EplHeadSubstepEnum);
739@@ -383,7 +385,7 @@
740 dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum, time); _assert_(dummy_input);
741 }
742 else{
743- dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time); _assert_(dummy_input);
744+ dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time,smb_averaging); _assert_(dummy_input);
745 }
746 surface_runoff_input=xDynamicCast<Input2*>(dummy_input); _assert_(surface_runoff_input);
747 }
Note: See TracBrowser for help on using the repository browser.