1 | /*!\file TransientInput.c
|
---|
2 | * \brief: implementation of the TransientInput object
|
---|
3 | */
|
---|
4 | /*Headers{{{*/
|
---|
5 | #ifdef HAVE_CONFIG_H
|
---|
6 | #include <config.h>
|
---|
7 | #else
|
---|
8 | #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
|
---|
9 | #endif
|
---|
10 |
|
---|
11 | #include <stdio.h>
|
---|
12 | #include <string.h>
|
---|
13 | #include "../objects.h"
|
---|
14 | #include "../../EnumDefinitions/EnumDefinitions.h"
|
---|
15 | #include "../../shared/shared.h"
|
---|
16 | #include "../../Container/Container.h"
|
---|
17 | #include "../../include/include.h"
|
---|
18 | /*}}}*/
|
---|
19 |
|
---|
20 | /*TransientInput constructors and destructor*/
|
---|
21 | /*FUNCTION TransientInput::TransientInput(){{{*/
|
---|
22 | TransientInput::TransientInput(){
|
---|
23 |
|
---|
24 | enum_type=UNDEF;
|
---|
25 | inputs=NULL;
|
---|
26 | this->numtimesteps=0;
|
---|
27 | this->parameters=NULL;
|
---|
28 | this->timesteps=NULL;
|
---|
29 |
|
---|
30 | }
|
---|
31 | /*}}}*/
|
---|
32 | /*FUNCTION TransientInput::TransientInput(int in_enum_type){{{*/
|
---|
33 | TransientInput::TransientInput(int in_enum_type)
|
---|
34 | {
|
---|
35 | /*Set Enum*/
|
---|
36 | enum_type=in_enum_type;
|
---|
37 |
|
---|
38 | /*Allocate values and timesteps, and copy: */
|
---|
39 | this->numtimesteps=0;
|
---|
40 | this->timesteps=NULL;
|
---|
41 | inputs = new Inputs();
|
---|
42 | this->parameters=NULL;
|
---|
43 |
|
---|
44 | }
|
---|
45 | /*}}}*/
|
---|
46 | /*FUNCTION TransientInput::~TransientInput{{{*/
|
---|
47 | TransientInput::~TransientInput(){
|
---|
48 | xDelete<IssmDouble>(this->timesteps);
|
---|
49 | this->timesteps=NULL;
|
---|
50 | this->numtimesteps=0;
|
---|
51 | parameters=NULL;
|
---|
52 | delete this->inputs;
|
---|
53 | return;
|
---|
54 | }
|
---|
55 | /*}}}*/
|
---|
56 |
|
---|
57 | /*Object virtual functions definitions:*/
|
---|
58 | /*FUNCTION TransientInput::Echo {{{*/
|
---|
59 | void TransientInput::Echo(void){
|
---|
60 | this->DeepEcho();
|
---|
61 | }
|
---|
62 | /*}}}*/
|
---|
63 | /*FUNCTION TransientInput::DeepEcho{{{*/
|
---|
64 | void TransientInput::DeepEcho(void){
|
---|
65 |
|
---|
66 | int i;
|
---|
67 |
|
---|
68 | _printLine_("TransientInput:");
|
---|
69 | _printLine_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")");
|
---|
70 | _printLine_(" numtimesteps: " << this->numtimesteps);
|
---|
71 | _printLine_("---inputs: ");
|
---|
72 | for(i=0;i<this->numtimesteps;i++){
|
---|
73 | _printLine_(" time: " << this->timesteps[i] << " ");
|
---|
74 | ((Input*)this->inputs->GetObjectByOffset(i))->Echo();
|
---|
75 | }
|
---|
76 | }
|
---|
77 | /*}}}*/
|
---|
78 | /*FUNCTION TransientInput::Id{{{*/
|
---|
79 | int TransientInput::Id(void){ return -1; }
|
---|
80 | /*}}}*/
|
---|
81 | /*FUNCTION TransientInput::MyRank{{{*/
|
---|
82 | int TransientInput::MyRank(void){
|
---|
83 | extern int my_rank;
|
---|
84 | return my_rank;
|
---|
85 | }
|
---|
86 | /*}}}*/
|
---|
87 | /*FUNCTION TransientInput::ObjectEnum{{{*/
|
---|
88 | int TransientInput::ObjectEnum(void){
|
---|
89 |
|
---|
90 | return TransientInputEnum;
|
---|
91 |
|
---|
92 | }
|
---|
93 | /*}}}*/
|
---|
94 | /*FUNCTION TransientInput::copy{{{*/
|
---|
95 | Object* TransientInput::copy() {
|
---|
96 |
|
---|
97 | TransientInput* output=NULL;
|
---|
98 |
|
---|
99 | output = new TransientInput();
|
---|
100 | output->enum_type=this->enum_type;
|
---|
101 | output->numtimesteps=this->numtimesteps;
|
---|
102 | output->timesteps=xNew<IssmDouble>(this->numtimesteps);
|
---|
103 | memcpy(output->timesteps,this->timesteps,this->numtimesteps*sizeof(IssmDouble));
|
---|
104 | output->inputs=(Inputs*)this->inputs->Copy();
|
---|
105 | output->parameters=this->parameters;
|
---|
106 |
|
---|
107 | return output;
|
---|
108 |
|
---|
109 | }
|
---|
110 | /*}}}*/
|
---|
111 |
|
---|
112 | /*TransientInput management*/
|
---|
113 | /*FUNCTION TransientInput::InstanceEnum{{{*/
|
---|
114 | int TransientInput::InstanceEnum(void){
|
---|
115 |
|
---|
116 | return this->enum_type;
|
---|
117 |
|
---|
118 | }
|
---|
119 | /*}}}*/
|
---|
120 | /*FUNCTION TransientInput::SpawnTriaInput{{{*/
|
---|
121 | Input* TransientInput::SpawnTriaInput(int* indices){
|
---|
122 |
|
---|
123 | /*output*/
|
---|
124 | TransientInput* outinput=NULL;
|
---|
125 |
|
---|
126 | /*Create new Transientinput (copy of current input)*/
|
---|
127 | outinput=new TransientInput();
|
---|
128 | outinput->enum_type=this->enum_type;
|
---|
129 | outinput->numtimesteps=this->numtimesteps;
|
---|
130 | outinput->timesteps=xNew<IssmDouble>(this->numtimesteps);
|
---|
131 | memcpy(outinput->timesteps,this->timesteps,this->numtimesteps*sizeof(IssmDouble));
|
---|
132 | outinput->inputs=(Inputs*)this->inputs->SpawnTriaInputs(indices);
|
---|
133 | outinput->parameters=this->parameters;
|
---|
134 |
|
---|
135 | /*Assign output*/
|
---|
136 | return outinput;
|
---|
137 |
|
---|
138 | }
|
---|
139 | /*}}}*/
|
---|
140 | /*FUNCTION TransientInput::SpawnResult{{{*/
|
---|
141 | ElementResult* TransientInput::SpawnResult(int step, IssmDouble time){
|
---|
142 |
|
---|
143 | ElementResult* elementresult=NULL;
|
---|
144 |
|
---|
145 | /*Ok, we want to spawn an ElementResult. We have the time, just get
|
---|
146 | *the correct values: */
|
---|
147 | Input* input=GetTimeInput(time);
|
---|
148 |
|
---|
149 | elementresult=input->SpawnResult(step,time);
|
---|
150 |
|
---|
151 | delete input;
|
---|
152 |
|
---|
153 | return elementresult;
|
---|
154 | }
|
---|
155 | /*}}}*/
|
---|
156 |
|
---|
157 | /*Object functions*/
|
---|
158 | /*FUNCTION TransientInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss){{{*/
|
---|
159 | void TransientInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss){
|
---|
160 | IssmDouble time;
|
---|
161 |
|
---|
162 | /*First, recover current time from parameters: */
|
---|
163 | this->parameters->FindParam(&time,TimeEnum);
|
---|
164 |
|
---|
165 | /*Retrieve interpolated values for this time step: */
|
---|
166 | Input* input=GetTimeInput(time);
|
---|
167 |
|
---|
168 | /*Call input function*/
|
---|
169 | input->GetInputValue(pvalue,gauss);
|
---|
170 |
|
---|
171 | delete input;
|
---|
172 | }
|
---|
173 | /*}}}*/
|
---|
174 | /*FUNCTION TransientInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){{{*/
|
---|
175 | void TransientInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){
|
---|
176 |
|
---|
177 | /*Retrieve interpolated values for this time step: */
|
---|
178 | Input* input=GetTimeInput(time);
|
---|
179 |
|
---|
180 | /*Call input function*/
|
---|
181 | input->GetInputValue(pvalue,gauss);
|
---|
182 |
|
---|
183 | delete input;
|
---|
184 | }
|
---|
185 | /*}}}*/
|
---|
186 | /*FUNCTION TransientInput::GetInputDerivativeValue(IssmDouble* p, IssmPDouble* xyz_list, GaussTria* gauss){{{*/
|
---|
187 | void TransientInput::GetInputDerivativeValue(IssmDouble* p, IssmPDouble* xyz_list, GaussTria* gauss){
|
---|
188 |
|
---|
189 | IssmDouble time;
|
---|
190 |
|
---|
191 | /*First, recover current time from parameters: */
|
---|
192 | parameters->FindParam(&time,TimeEnum);
|
---|
193 |
|
---|
194 | /*Retrieve interpolated values for this time step: */
|
---|
195 | Input* input=GetTimeInput(time);
|
---|
196 |
|
---|
197 | /*Call input function*/
|
---|
198 | input->GetInputDerivativeValue(p,xyz_list,gauss);
|
---|
199 |
|
---|
200 | delete input;
|
---|
201 |
|
---|
202 | }
|
---|
203 | /*}}}*/
|
---|
204 | /*FUNCTION TransientInput::ChangeEnum{{{*/
|
---|
205 | void TransientInput::ChangeEnum(int newenumtype){
|
---|
206 | this->enum_type=newenumtype;
|
---|
207 | }
|
---|
208 | /*}}}*/
|
---|
209 | /*FUNCTION TransientInput::GetInputAverage{{{*/
|
---|
210 | void TransientInput::GetInputAverage(IssmDouble* pvalue){
|
---|
211 |
|
---|
212 | IssmDouble time;
|
---|
213 |
|
---|
214 | /*First, recover current time from parameters: */
|
---|
215 | parameters->FindParam(&time,TimeEnum);
|
---|
216 |
|
---|
217 | /*Retrieve interpolated values for this time step: */
|
---|
218 | Input* input=GetTimeInput(time);
|
---|
219 |
|
---|
220 | /*Call input function*/
|
---|
221 | input->GetInputAverage(pvalue);
|
---|
222 |
|
---|
223 | delete input;
|
---|
224 |
|
---|
225 | }
|
---|
226 | /*}}}*/
|
---|
227 |
|
---|
228 | /*Intermediary*/
|
---|
229 | /*FUNCTION TransientInput::AddTimeInput{{{*/
|
---|
230 | void TransientInput::AddTimeInput(Input* input,IssmDouble time){
|
---|
231 |
|
---|
232 | /*insert values at time step: */
|
---|
233 | if (this->numtimesteps>0 && time<=this->timesteps[this->numtimesteps-1]) _assert_("timestep values must increase sequentially");
|
---|
234 |
|
---|
235 | //copy timesteps, add the new time, delete previous timesteps, and add the new input: inputs->AddObject(input);
|
---|
236 | IssmDouble* old_timesteps=NULL;
|
---|
237 |
|
---|
238 | if (this->numtimesteps > 0){
|
---|
239 | old_timesteps=xNew<IssmDouble>(this->numtimesteps);
|
---|
240 | memcpy(old_timesteps,this->timesteps,this->numtimesteps*sizeof(IssmDouble));
|
---|
241 | xDelete<IssmDouble>(this->timesteps);
|
---|
242 | }
|
---|
243 |
|
---|
244 | this->numtimesteps=this->numtimesteps+1;
|
---|
245 | this->timesteps=xNew<IssmDouble>(this->numtimesteps);
|
---|
246 |
|
---|
247 | if (this->numtimesteps > 1){
|
---|
248 | memcpy(this->timesteps,old_timesteps,(this->numtimesteps-1)*sizeof(IssmDouble));
|
---|
249 | xDelete<IssmDouble>(old_timesteps);
|
---|
250 | }
|
---|
251 |
|
---|
252 | /*go ahead and plug: */
|
---|
253 | this->timesteps[this->numtimesteps-1]=time;
|
---|
254 | inputs->AddObject(input);
|
---|
255 |
|
---|
256 | }
|
---|
257 | /*}}}*/
|
---|
258 | /*FUNCTION TransientInput::SquareMin{{{*/
|
---|
259 | void TransientInput::SquareMin(IssmDouble* psquaremin, bool process_units,Parameters* parameters){
|
---|
260 |
|
---|
261 | IssmDouble time;
|
---|
262 |
|
---|
263 | /*First, recover current time from parameters: */
|
---|
264 | parameters->FindParam(&time,TimeEnum);
|
---|
265 |
|
---|
266 | /*Retrieve interpolated values for this time step: */
|
---|
267 | Input* input=GetTimeInput(time);
|
---|
268 |
|
---|
269 | /*Call input function*/
|
---|
270 | input->SquareMin(psquaremin,process_units,parameters);
|
---|
271 |
|
---|
272 | delete input;
|
---|
273 |
|
---|
274 | }
|
---|
275 | /*}}}*/
|
---|
276 | /*FUNCTION TransientInput::InfinityNorm{{{*/
|
---|
277 | IssmDouble TransientInput::InfinityNorm(void){
|
---|
278 |
|
---|
279 | IssmDouble time;
|
---|
280 | IssmDouble infnorm;
|
---|
281 |
|
---|
282 | /*First, recover current time from parameters: */
|
---|
283 | parameters->FindParam(&time,TimeEnum);
|
---|
284 |
|
---|
285 | /*Retrieve interpolated values for this time step: */
|
---|
286 | Input* input=GetTimeInput(time);
|
---|
287 |
|
---|
288 | /*Call input function*/
|
---|
289 | infnorm=input->InfinityNorm();
|
---|
290 |
|
---|
291 | /*Clean-up and return*/
|
---|
292 | delete input;
|
---|
293 | return infnorm;
|
---|
294 | }
|
---|
295 | /*}}}*/
|
---|
296 | /*FUNCTION TransientInput::Max{{{*/
|
---|
297 | IssmDouble TransientInput::Max(void){
|
---|
298 |
|
---|
299 | IssmDouble time;
|
---|
300 | IssmDouble max;
|
---|
301 |
|
---|
302 | /*First, recover current time from parameters: */
|
---|
303 | parameters->FindParam(&time,TimeEnum);
|
---|
304 |
|
---|
305 | /*Retrieve interpolated values for this time step: */
|
---|
306 | Input* input=GetTimeInput(time);
|
---|
307 |
|
---|
308 | /*Call input function*/
|
---|
309 | max=input->Max();
|
---|
310 |
|
---|
311 | delete input;
|
---|
312 |
|
---|
313 | return max;
|
---|
314 | }
|
---|
315 | /*}}}*/
|
---|
316 | /*FUNCTION TransientInput::MaxAbs{{{*/
|
---|
317 | IssmDouble TransientInput::MaxAbs(void){
|
---|
318 |
|
---|
319 | IssmDouble time;
|
---|
320 | IssmDouble maxabs;
|
---|
321 |
|
---|
322 | /*First, recover current time from parameters: */
|
---|
323 | parameters->FindParam(&time,TimeEnum);
|
---|
324 |
|
---|
325 | /*Retrieve interpolated values for this time step: */
|
---|
326 | Input* input=GetTimeInput(time);
|
---|
327 |
|
---|
328 | /*Call input function*/
|
---|
329 | maxabs=input->MaxAbs();
|
---|
330 |
|
---|
331 | /*Clean-up and return*/
|
---|
332 | delete input;
|
---|
333 | return maxabs;
|
---|
334 |
|
---|
335 | }
|
---|
336 | /*}}}*/
|
---|
337 | /*FUNCTION TransientInput::Min{{{*/
|
---|
338 | IssmDouble TransientInput::Min(void){
|
---|
339 |
|
---|
340 | IssmDouble time;
|
---|
341 | IssmDouble min;
|
---|
342 |
|
---|
343 | /*First, recover current time from parameters: */
|
---|
344 | parameters->FindParam(&time,TimeEnum);
|
---|
345 |
|
---|
346 | /*Retrieve interpolated values for this time step: */
|
---|
347 | Input* input=GetTimeInput(time);
|
---|
348 |
|
---|
349 | /*Call input function*/
|
---|
350 | min=input->Min();
|
---|
351 |
|
---|
352 | /*Clean-up and return*/
|
---|
353 | delete input;
|
---|
354 | return min;
|
---|
355 |
|
---|
356 | }
|
---|
357 | /*}}}*/
|
---|
358 | /*FUNCTION TransientInput::MinAbs{{{*/
|
---|
359 | IssmDouble TransientInput::MinAbs(void){
|
---|
360 |
|
---|
361 | IssmDouble time;
|
---|
362 | IssmDouble minabs;
|
---|
363 |
|
---|
364 | /*First, recover current time from parameters: */
|
---|
365 | parameters->FindParam(&time,TimeEnum);
|
---|
366 |
|
---|
367 | /*Retrieve interpolated values for this time step: */
|
---|
368 | Input* input=GetTimeInput(time);
|
---|
369 |
|
---|
370 | /*Call input function*/
|
---|
371 | minabs=input->MinAbs();
|
---|
372 |
|
---|
373 | /*Clean-up and return*/
|
---|
374 | delete input;
|
---|
375 | return minabs;
|
---|
376 | }
|
---|
377 | /*}}}*/
|
---|
378 | /*FUNCTION TransientInput::GetVectorFromInputs{{{*/
|
---|
379 | void TransientInput::GetVectorFromInputs(Vector* vector,int* doflist){
|
---|
380 |
|
---|
381 | IssmDouble time;
|
---|
382 |
|
---|
383 | /*First, recover current time from parameters: */
|
---|
384 | parameters->FindParam(&time,TimeEnum);
|
---|
385 |
|
---|
386 | /*Retrieve interpolated values for this time step: */
|
---|
387 | Input* input=GetTimeInput(time);
|
---|
388 |
|
---|
389 | /*Call input function*/
|
---|
390 | input->GetVectorFromInputs(vector,doflist);
|
---|
391 |
|
---|
392 | delete input;
|
---|
393 |
|
---|
394 | } /*}}}*/
|
---|
395 | /*FUNCTION TransientInput::GetTimeInput{{{*/
|
---|
396 | Input* TransientInput::GetTimeInput(IssmDouble intime){
|
---|
397 |
|
---|
398 | int i,j;
|
---|
399 | IssmDouble deltat;
|
---|
400 | IssmDouble alpha1,alpha2;
|
---|
401 | bool found=false;
|
---|
402 | Input* input=NULL;
|
---|
403 | Input* input1=NULL;
|
---|
404 | Input* input2=NULL;
|
---|
405 |
|
---|
406 | /*Ok, we have the time, go through the timesteps, and figure out which interval we
|
---|
407 | *fall within. Then interpolate the values on this interval: */
|
---|
408 | if(intime<this->timesteps[0]){
|
---|
409 | /*get values for the first time: */
|
---|
410 | input=(Input*)((Input*)this->inputs->GetObjectByOffset(0))->copy();
|
---|
411 | found=true;
|
---|
412 | }
|
---|
413 | else if(intime>this->timesteps[this->numtimesteps-1]){
|
---|
414 | /*get values for the last time: */
|
---|
415 | input=(Input*)((Input*)this->inputs->GetObjectByOffset(numtimesteps-1))->copy();
|
---|
416 | found=true;
|
---|
417 | }
|
---|
418 | else{
|
---|
419 | /*Find which interval we fall within: */
|
---|
420 | for(i=0;i<this->numtimesteps;i++){
|
---|
421 | if(intime==this->timesteps[i]){
|
---|
422 | /*We are right on one step time: */
|
---|
423 | input=(Input*)((Input*)this->inputs->GetObjectByOffset(i))->copy();
|
---|
424 | found=true;
|
---|
425 | break; //we are done with the time interpolation.
|
---|
426 | }
|
---|
427 | else{
|
---|
428 | if(this->timesteps[i]<intime && intime<this->timesteps[i+1]){
|
---|
429 | /*ok, we have the interval ]i:i+1[. Interpolate linearly for now: */
|
---|
430 | deltat=this->timesteps[i+1]-this->timesteps[i];
|
---|
431 | alpha2=(intime-this->timesteps[i])/deltat;
|
---|
432 | alpha1=(1-alpha2);
|
---|
433 |
|
---|
434 | input1=(Input*)this->inputs->GetObjectByOffset(i);
|
---|
435 | input2=(Input*)this->inputs->GetObjectByOffset(i+1);
|
---|
436 |
|
---|
437 | input=(Input*)input1->copy();
|
---|
438 | input->Scale(alpha1);
|
---|
439 | input->AXPY(input2,alpha2);
|
---|
440 |
|
---|
441 | found=true;
|
---|
442 | break;
|
---|
443 | }
|
---|
444 | else continue; //keep looking on the next interval
|
---|
445 | }
|
---|
446 | }
|
---|
447 | }
|
---|
448 | if(!found)_error2_("did not find time interval on which to interpolate forcing values!");
|
---|
449 |
|
---|
450 | /*Assign output pointer*/
|
---|
451 | return input;
|
---|
452 | }
|
---|
453 | /*}}}*/
|
---|
454 | /*FUNCTION TransientInput::Configure{{{*/
|
---|
455 | void TransientInput::Configure(Parameters* parameters){
|
---|
456 | this->parameters=parameters;
|
---|
457 | }
|
---|
458 | /*}}}*/
|
---|