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

/*!\file TransientInput.c
 * \brief: implementation of the TransientInput object
 */
/*Headers{{{*/
#ifdef HAVE_CONFIG_H
        #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <stdio.h>
#include <string.h>
#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<double>(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;

        printf("TransientInput:\n");
        printf("   enum: %i (%s)\n",this->enum_type,EnumToStringx(this->enum_type));
        printf("   numtimesteps: %i\n",this->numtimesteps);
        printf("---inputs: \n"); 
        for(i=0;i<this->numtimesteps;i++){
                printf("   time: %g  \n",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<IssmPDouble>(this->numtimesteps);
   memcpy(output->timesteps,this->timesteps,this->numtimesteps*sizeof(IssmPDouble));
        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<IssmPDouble>(this->numtimesteps);
        memcpy(outinput->timesteps,this->timesteps,this->numtimesteps*sizeof(IssmPDouble));
        outinput->inputs=(Inputs*)this->inputs->SpawnTriaInputs(indices);
        outinput->parameters=this->parameters;

        /*Assign output*/
        return outinput;

}
/*}}}*/
/*FUNCTION TransientInput::SpawnResult{{{*/
ElementResult* TransientInput::SpawnResult(int step, IssmPDouble 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(IssmPDouble* pvalue,GaussTria* gauss){{{*/
void TransientInput::GetInputValue(IssmPDouble* pvalue,GaussTria* gauss){
        IssmPDouble 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(IssmPDouble* pvalue,GaussTria* gauss,IssmPDouble time){{{*/
void TransientInput::GetInputValue(IssmPDouble* pvalue,GaussTria* gauss,IssmPDouble time){

        /*Retrieve interpolated values for this time step: */
        Input* input=GetTimeInput(time);

        /*Call input function*/
        input->GetInputValue(pvalue,gauss);

        delete input;
}
/*}}}*/
/*FUNCTION TransientInput::GetInputDerivativeValue(IssmPDouble* p, IssmPDouble* xyz_list, GaussTria* gauss){{{*/
void TransientInput::GetInputDerivativeValue(IssmPDouble* p, IssmPDouble* xyz_list, GaussTria* gauss){

        IssmPDouble 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(IssmPDouble* pvalue){
        
        IssmPDouble 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,IssmPDouble 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);
        IssmPDouble* old_timesteps=NULL;

        if (this->numtimesteps > 0){
                old_timesteps=xNew<IssmPDouble>(this->numtimesteps);
                memcpy(old_timesteps,this->timesteps,this->numtimesteps*sizeof(IssmPDouble));
                xDelete<double>(this->timesteps); 
        }

        this->numtimesteps=this->numtimesteps+1;
        this->timesteps=xNew<IssmPDouble>(this->numtimesteps);

        if (this->numtimesteps > 1){
                memcpy(this->timesteps,old_timesteps,(this->numtimesteps-1)*sizeof(IssmPDouble));
                xDelete<double>(old_timesteps);
        }

        /*go ahead and plug: */
        this->timesteps[this->numtimesteps-1]=time;
        inputs->AddObject(input);

}
/*}}}*/
/*FUNCTION TransientInput::SquareMin{{{*/
void TransientInput::SquareMin(IssmPDouble* psquaremin, bool process_units,Parameters* parameters){

        IssmPDouble 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{{{*/
IssmPDouble TransientInput::InfinityNorm(void){

        IssmPDouble time;
        IssmPDouble 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{{{*/
IssmPDouble TransientInput::Max(void){

        IssmPDouble time;
        IssmPDouble 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{{{*/
IssmPDouble TransientInput::MaxAbs(void){

        IssmPDouble time;
        IssmPDouble 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{{{*/
IssmPDouble TransientInput::Min(void){

        IssmPDouble time;
        IssmPDouble 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{{{*/
IssmPDouble TransientInput::MinAbs(void){

        IssmPDouble time;
        IssmPDouble 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){

        IssmPDouble 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(IssmPDouble intime){

        int     i,j;
        IssmPDouble  deltat;
        IssmPDouble  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(intime<this->timesteps[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;i<this->numtimesteps;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]<intime && intime<this->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;
}
/*}}}*/