source: issm/trunk-jpl/src/c/objects/Inputs/TransientInput.cpp@ 12530

Last change on this file since 12530 was 12530, checked in by utke, 13 years ago

in the first round the assumption was made that the inputs remain passive but in the second round after discussion with M. Morlighem all inputs were consistently made active with the exception of the xyz_list arguments

File size: 11.3 KB
Line 
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(){{{*/
22TransientInput::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){{{*/
33TransientInput::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{{{*/
47TransientInput::~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 {{{*/
59void TransientInput::Echo(void){
60 this->DeepEcho();
61}
62/*}}}*/
63/*FUNCTION TransientInput::DeepEcho{{{*/
64void 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{{{*/
79int TransientInput::Id(void){ return -1; }
80/*}}}*/
81/*FUNCTION TransientInput::MyRank{{{*/
82int TransientInput::MyRank(void){
83 extern int my_rank;
84 return my_rank;
85}
86/*}}}*/
87/*FUNCTION TransientInput::ObjectEnum{{{*/
88int TransientInput::ObjectEnum(void){
89
90 return TransientInputEnum;
91
92}
93/*}}}*/
94/*FUNCTION TransientInput::copy{{{*/
95Object* 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{{{*/
114int TransientInput::InstanceEnum(void){
115
116 return this->enum_type;
117
118}
119/*}}}*/
120/*FUNCTION TransientInput::SpawnTriaInput{{{*/
121Input* 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{{{*/
141ElementResult* 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){{{*/
159void 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){{{*/
175void 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){{{*/
187void 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{{{*/
205void TransientInput::ChangeEnum(int newenumtype){
206 this->enum_type=newenumtype;
207}
208/*}}}*/
209/*FUNCTION TransientInput::GetInputAverage{{{*/
210void 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{{{*/
230void 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{{{*/
259void 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{{{*/
277IssmDouble 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{{{*/
297IssmDouble 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{{{*/
317IssmDouble 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{{{*/
338IssmDouble 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{{{*/
359IssmDouble 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{{{*/
379void 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{{{*/
396Input* 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{{{*/
455void TransientInput::Configure(Parameters* parameters){
456 this->parameters=parameters;
457}
458/*}}}*/
Note: See TracBrowser for help on using the repository browser.