/*!\file TransientInput.c * \brief: implementation of the TransientInput object */ /*Headers{{{*/ #ifdef HAVE_CONFIG_H #include #else #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" #endif #include #include #include "../objects.h" #include "../../EnumDefinitions/EnumDefinitions.h" #include "../../shared/shared.h" #include "../../Container/Container.h" #include "../../include/include.h" /*}}}*/ /*TransientInput constructors and destructor*/ /*FUNCTION TransientInput::TransientInput(){{{*/ TransientInput::TransientInput(){ enum_type=UNDEF; inputs=NULL; this->numtimesteps=0; this->parameters=NULL; this->timesteps=NULL; } /*}}}*/ /*FUNCTION TransientInput::TransientInput(int in_enum_type){{{*/ TransientInput::TransientInput(int in_enum_type) { /*Set Enum*/ enum_type=in_enum_type; /*Allocate values and timesteps, and copy: */ this->numtimesteps=0; this->timesteps=NULL; inputs = new Inputs(); this->parameters=NULL; } /*}}}*/ /*FUNCTION TransientInput::~TransientInput{{{*/ TransientInput::~TransientInput(){ xDelete(this->timesteps); this->timesteps=NULL; this->numtimesteps=0; parameters=NULL; delete this->inputs; return; } /*}}}*/ /*Object virtual functions definitions:*/ /*FUNCTION TransientInput::Echo {{{*/ void TransientInput::Echo(void){ this->DeepEcho(); } /*}}}*/ /*FUNCTION TransientInput::DeepEcho{{{*/ void TransientInput::DeepEcho(void){ int i; _printLine_("TransientInput:"); _printLine_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")"); _printLine_(" numtimesteps: " << this->numtimesteps); _printLine_("---inputs: "); for(i=0;inumtimesteps;i++){ _printLine_(" time: " << this->timesteps[i] << " "); ((Input*)this->inputs->GetObjectByOffset(i))->Echo(); } } /*}}}*/ /*FUNCTION TransientInput::Id{{{*/ int TransientInput::Id(void){ return -1; } /*}}}*/ /*FUNCTION TransientInput::MyRank{{{*/ int TransientInput::MyRank(void){ extern int my_rank; return my_rank; } /*}}}*/ /*FUNCTION TransientInput::ObjectEnum{{{*/ int TransientInput::ObjectEnum(void){ return TransientInputEnum; } /*}}}*/ /*FUNCTION TransientInput::copy{{{*/ Object* TransientInput::copy() { TransientInput* output=NULL; output = new TransientInput(); output->enum_type=this->enum_type; output->numtimesteps=this->numtimesteps; output->timesteps=xNew(this->numtimesteps); memcpy(output->timesteps,this->timesteps,this->numtimesteps*sizeof(IssmDouble)); output->inputs=(Inputs*)this->inputs->Copy(); output->parameters=this->parameters; return output; } /*}}}*/ /*TransientInput management*/ /*FUNCTION TransientInput::InstanceEnum{{{*/ int TransientInput::InstanceEnum(void){ return this->enum_type; } /*}}}*/ /*FUNCTION TransientInput::SpawnTriaInput{{{*/ Input* TransientInput::SpawnTriaInput(int* indices){ /*output*/ TransientInput* outinput=NULL; /*Create new Transientinput (copy of current input)*/ outinput=new TransientInput(); outinput->enum_type=this->enum_type; outinput->numtimesteps=this->numtimesteps; outinput->timesteps=xNew(this->numtimesteps); memcpy(outinput->timesteps,this->timesteps,this->numtimesteps*sizeof(IssmDouble)); outinput->inputs=(Inputs*)this->inputs->SpawnTriaInputs(indices); outinput->parameters=this->parameters; /*Assign output*/ return outinput; } /*}}}*/ /*FUNCTION TransientInput::SpawnResult{{{*/ ElementResult* TransientInput::SpawnResult(int step, IssmDouble time){ ElementResult* elementresult=NULL; /*Ok, we want to spawn an ElementResult. We have the time, just get *the correct values: */ Input* input=GetTimeInput(time); elementresult=input->SpawnResult(step,time); delete input; return elementresult; } /*}}}*/ /*Object functions*/ /*FUNCTION TransientInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss){{{*/ void TransientInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss){ IssmDouble time; /*First, recover current time from parameters: */ this->parameters->FindParam(&time,TimeEnum); /*Retrieve interpolated values for this time step: */ Input* input=GetTimeInput(time); /*Call input function*/ input->GetInputValue(pvalue,gauss); delete input; } /*}}}*/ /*FUNCTION TransientInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){{{*/ void TransientInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){ /*Retrieve interpolated values for this time step: */ Input* input=GetTimeInput(time); /*Call input function*/ input->GetInputValue(pvalue,gauss); delete input; } /*}}}*/ /*FUNCTION TransientInput::GetInputDerivativeValue(IssmDouble* p, IssmPDouble* xyz_list, GaussTria* gauss){{{*/ void TransientInput::GetInputDerivativeValue(IssmDouble* p, IssmPDouble* xyz_list, GaussTria* gauss){ IssmDouble time; /*First, recover current time from parameters: */ parameters->FindParam(&time,TimeEnum); /*Retrieve interpolated values for this time step: */ Input* input=GetTimeInput(time); /*Call input function*/ input->GetInputDerivativeValue(p,xyz_list,gauss); delete input; } /*}}}*/ /*FUNCTION TransientInput::ChangeEnum{{{*/ void TransientInput::ChangeEnum(int newenumtype){ this->enum_type=newenumtype; } /*}}}*/ /*FUNCTION TransientInput::GetInputAverage{{{*/ void TransientInput::GetInputAverage(IssmDouble* pvalue){ IssmDouble time; /*First, recover current time from parameters: */ parameters->FindParam(&time,TimeEnum); /*Retrieve interpolated values for this time step: */ Input* input=GetTimeInput(time); /*Call input function*/ input->GetInputAverage(pvalue); delete input; } /*}}}*/ /*Intermediary*/ /*FUNCTION TransientInput::AddTimeInput{{{*/ void TransientInput::AddTimeInput(Input* input,IssmDouble time){ /*insert values at time step: */ if (this->numtimesteps>0 && time<=this->timesteps[this->numtimesteps-1]) _assert_("timestep values must increase sequentially"); //copy timesteps, add the new time, delete previous timesteps, and add the new input: inputs->AddObject(input); IssmDouble* old_timesteps=NULL; if (this->numtimesteps > 0){ old_timesteps=xNew(this->numtimesteps); memcpy(old_timesteps,this->timesteps,this->numtimesteps*sizeof(IssmDouble)); xDelete(this->timesteps); } this->numtimesteps=this->numtimesteps+1; this->timesteps=xNew(this->numtimesteps); if (this->numtimesteps > 1){ memcpy(this->timesteps,old_timesteps,(this->numtimesteps-1)*sizeof(IssmDouble)); xDelete(old_timesteps); } /*go ahead and plug: */ this->timesteps[this->numtimesteps-1]=time; inputs->AddObject(input); } /*}}}*/ /*FUNCTION TransientInput::SquareMin{{{*/ void TransientInput::SquareMin(IssmDouble* psquaremin, bool process_units,Parameters* parameters){ IssmDouble time; /*First, recover current time from parameters: */ parameters->FindParam(&time,TimeEnum); /*Retrieve interpolated values for this time step: */ Input* input=GetTimeInput(time); /*Call input function*/ input->SquareMin(psquaremin,process_units,parameters); delete input; } /*}}}*/ /*FUNCTION TransientInput::InfinityNorm{{{*/ IssmDouble TransientInput::InfinityNorm(void){ IssmDouble time; IssmDouble infnorm; /*First, recover current time from parameters: */ parameters->FindParam(&time,TimeEnum); /*Retrieve interpolated values for this time step: */ Input* input=GetTimeInput(time); /*Call input function*/ infnorm=input->InfinityNorm(); /*Clean-up and return*/ delete input; return infnorm; } /*}}}*/ /*FUNCTION TransientInput::Max{{{*/ IssmDouble TransientInput::Max(void){ IssmDouble time; IssmDouble max; /*First, recover current time from parameters: */ parameters->FindParam(&time,TimeEnum); /*Retrieve interpolated values for this time step: */ Input* input=GetTimeInput(time); /*Call input function*/ max=input->Max(); delete input; return max; } /*}}}*/ /*FUNCTION TransientInput::MaxAbs{{{*/ IssmDouble TransientInput::MaxAbs(void){ IssmDouble time; IssmDouble maxabs; /*First, recover current time from parameters: */ parameters->FindParam(&time,TimeEnum); /*Retrieve interpolated values for this time step: */ Input* input=GetTimeInput(time); /*Call input function*/ maxabs=input->MaxAbs(); /*Clean-up and return*/ delete input; return maxabs; } /*}}}*/ /*FUNCTION TransientInput::Min{{{*/ IssmDouble TransientInput::Min(void){ IssmDouble time; IssmDouble min; /*First, recover current time from parameters: */ parameters->FindParam(&time,TimeEnum); /*Retrieve interpolated values for this time step: */ Input* input=GetTimeInput(time); /*Call input function*/ min=input->Min(); /*Clean-up and return*/ delete input; return min; } /*}}}*/ /*FUNCTION TransientInput::MinAbs{{{*/ IssmDouble TransientInput::MinAbs(void){ IssmDouble time; IssmDouble minabs; /*First, recover current time from parameters: */ parameters->FindParam(&time,TimeEnum); /*Retrieve interpolated values for this time step: */ Input* input=GetTimeInput(time); /*Call input function*/ minabs=input->MinAbs(); /*Clean-up and return*/ delete input; return minabs; } /*}}}*/ /*FUNCTION TransientInput::GetVectorFromInputs{{{*/ void TransientInput::GetVectorFromInputs(Vector* vector,int* doflist){ IssmDouble time; /*First, recover current time from parameters: */ parameters->FindParam(&time,TimeEnum); /*Retrieve interpolated values for this time step: */ Input* input=GetTimeInput(time); /*Call input function*/ input->GetVectorFromInputs(vector,doflist); delete input; } /*}}}*/ /*FUNCTION TransientInput::GetTimeInput{{{*/ Input* TransientInput::GetTimeInput(IssmDouble intime){ int i,j; IssmDouble deltat; IssmDouble alpha1,alpha2; bool found=false; Input* input=NULL; Input* input1=NULL; Input* input2=NULL; /*Ok, we have the time, go through the timesteps, and figure out which interval we *fall within. Then interpolate the values on this interval: */ if(intimetimesteps[0]){ /*get values for the first time: */ input=(Input*)((Input*)this->inputs->GetObjectByOffset(0))->copy(); found=true; } else if(intime>this->timesteps[this->numtimesteps-1]){ /*get values for the last time: */ input=(Input*)((Input*)this->inputs->GetObjectByOffset(numtimesteps-1))->copy(); found=true; } else{ /*Find which interval we fall within: */ for(i=0;inumtimesteps;i++){ if(intime==this->timesteps[i]){ /*We are right on one step time: */ input=(Input*)((Input*)this->inputs->GetObjectByOffset(i))->copy(); found=true; break; //we are done with the time interpolation. } else{ if(this->timesteps[i]timesteps[i+1]){ /*ok, we have the interval ]i:i+1[. Interpolate linearly for now: */ deltat=this->timesteps[i+1]-this->timesteps[i]; alpha2=(intime-this->timesteps[i])/deltat; alpha1=(1-alpha2); input1=(Input*)this->inputs->GetObjectByOffset(i); input2=(Input*)this->inputs->GetObjectByOffset(i+1); input=(Input*)input1->copy(); input->Scale(alpha1); input->AXPY(input2,alpha2); found=true; break; } else continue; //keep looking on the next interval } } } if(!found)_error2_("did not find time interval on which to interpolate forcing values!"); /*Assign output pointer*/ return input; } /*}}}*/ /*FUNCTION TransientInput::Configure{{{*/ void TransientInput::Configure(Parameters* parameters){ this->parameters=parameters; } /*}}}*/