TransientInput2.cpp

Go to the documentation of this file.

 
/*Headers*/
#ifdef HAVE_CONFIG_H
   #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif
 
#include "./TransientInput2.h"
#include "./TriaInput2.h"
#include "./PentaInput2.h"
#include "../../shared/shared.h"
#include "../Params/Parameters.h"
 
/*TransientInput2 constructors and destructor*/
TransientInput2::TransientInput2(){/*{{{*/
 
   enum_type=UNDEF;
   inputs=NULL;
   this->numtimesteps=0;
   this->parameters=NULL;
   this->timesteps=NULL;
 
   this->current_input=NULL;
   this->current_step=-1;
 
}
/*}}}*/
TransientInput2::TransientInput2(int in_enum_type,int nbe,int nbv,IssmDouble* timesin,int N){/*{{{*/
 
   /*Set Enum*/
   this->enum_type=in_enum_type;
   this->numberofelements_local = nbe;
   this->numberofvertices_local = nbv;
 
   /*Allocate values and timesteps, and copy: */
   _assert_(N>=0 && N<1e6);
   this->numtimesteps=N;
   if(N>0){
      this->timesteps=xNew<IssmDouble>(N);
      xMemCpy(this->timesteps,timesin,N);
 
      this->inputs     = xNew<Input2*>(N);
      for(int i=0;i<N;i++) this->inputs[i] = NULL;
   }
   else{
      this->timesteps=0;
      this->inputs   =0;
   }
   this->parameters = NULL;
   this->current_input=NULL;
   this->current_step=-1;
}
/*}}}*/
TransientInput2::~TransientInput2(){/*{{{*/
 
   for(int i=0;i<this->numtimesteps;i++){
      delete this->inputs[i];
   }
   xDelete<Input2*>(this->inputs);
   xDelete<IssmDouble>(this->timesteps);
 
   if(this->current_input) delete this->current_input;
}
/*}}}*/
 
/*Object virtual functions definitions:*/
Input2* TransientInput2::copy() {/*{{{*/
 
   TransientInput2* output=NULL;
 
   output = new TransientInput2();
   output->enum_type=this->enum_type;
   output->numtimesteps=this->numtimesteps;
   if(this->numtimesteps>0){
      output->timesteps=xNew<IssmDouble>(this->numtimesteps);
      xMemCpy(output->timesteps,this->timesteps,this->numtimesteps);
      output->inputs = xNew<Input2*>(this->numtimesteps);
      for(int i=0;i<this->numtimesteps;i++){
         output->inputs[i] = this->inputs[i]->copy();
      }
   }
   output->parameters=this->parameters;
 
   return output;
}/*}}}*/
void TransientInput2::DeepEcho(void){/*{{{*/
 
   int i;
 
   _printf_("TransientInput2:\n");
   _printf_("   enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
   _printf_("   numtimesteps: " << this->numtimesteps << "\n");
   _printf_("---inputs: \n");
   for(i=0;i<this->numtimesteps;i++){
      _printf_("   time: " << this->timesteps[i]<<"  ");
      if(this->inputs[i]) this->inputs[i]->Echo();
      else                _printf_(" NOT SET! \n");
   }
}
/*}}}*/
void TransientInput2::Configure(Parameters* params){/*{{{*/
   this->parameters=params;
}
/*}}}*/
void TransientInput2::Echo(void){/*{{{*/
   this->DeepEcho();
}
/*}}}*/
int  TransientInput2::Id(void){ return -1; }/*{{{*/
/*}}}*/
void TransientInput2::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
 
   if (marshall_direction == MARSHALLING_BACKWARD){
      _error_("not implmented");
      //inputs = new Inputs();
   }
 
   MARSHALLING_ENUM(TransientInput2Enum);
 
   MARSHALLING(enum_type);
   MARSHALLING(numtimesteps);
   MARSHALLING_DYNAMIC(this->timesteps,IssmDouble,numtimesteps);
   //inputs->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
   _error_("not implemented");
}
/*}}}*/
int  TransientInput2::ObjectEnum(void){/*{{{*/
 
   return TransientInput2Enum;
 
}
/*}}}*/
 
/*Intermediary*/
void TransientInput2::AddTriaTimeInput(IssmDouble time,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/
 
   /*Check whether this is the last time step that we have*/
   if(this->numtimesteps){
      if(fabs(this->timesteps[this->numtimesteps-1]-time)<1.0e-5){
         this->AddTriaTimeInput(this->numtimesteps-1,numindices,indices,values_in,interp_in);
         return;
      }
   }
 
   /*This is a new time step! we need to add it to the list*/
   if(this->numtimesteps>0 && time<this->timesteps[this->numtimesteps-1]) _error_("timestep values must increase sequentially");
 
   IssmDouble *old_timesteps = NULL;
   Input2    **old_inputs    = NULL;
   if (this->numtimesteps > 0){
      old_timesteps=xNew<IssmDouble>(this->numtimesteps);
      xMemCpy(old_timesteps,this->timesteps,this->numtimesteps);
      xDelete<IssmDouble>(this->timesteps);
      old_inputs=xNew<Input2*>(this->numtimesteps);
      xMemCpy(old_inputs,this->inputs,this->numtimesteps);
      xDelete<Input2*>(this->inputs);
   }
 
   this->numtimesteps=this->numtimesteps+1;
   this->timesteps=xNew<IssmDouble>(this->numtimesteps);
   this->inputs   = xNew<Input2*>(this->numtimesteps);
 
   if (this->numtimesteps > 1){
      xMemCpy(this->inputs,old_inputs,this->numtimesteps-1);
      xMemCpy(this->timesteps,old_timesteps,this->numtimesteps-1);
      xDelete(old_timesteps);
      xDelete<Input2*>(old_inputs);
   }
 
   /*go ahead and plug: */
   this->timesteps[this->numtimesteps-1] = time;
   this->inputs[this->numtimesteps-1]    = NULL;
   this->AddTriaTimeInput(this->numtimesteps-1,numindices,indices,values_in,interp_in);
 
}
/*}}}*/
void TransientInput2::AddPentaTimeInput(IssmDouble time,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/
 
   /*Check whether this is the last time step that we have*/
   if(this->numtimesteps){
      if(fabs(this->timesteps[this->numtimesteps-1]-time)<1.0e-5){
         this->AddPentaTimeInput(this->numtimesteps-1,numindices,indices,values_in,interp_in);
         return;
      }
   }
 
   /*This is a new time step! we need to add it to the list*/
   if(this->numtimesteps>0 && time<this->timesteps[this->numtimesteps-1]) _error_("timestep values must increase sequentially");
 
   IssmDouble *old_timesteps = NULL;
   Input2    **old_inputs    = NULL;
   if (this->numtimesteps > 0){
      old_timesteps=xNew<IssmDouble>(this->numtimesteps);
      xMemCpy(old_timesteps,this->timesteps,this->numtimesteps);
      xDelete<IssmDouble>(this->timesteps);
      old_inputs=xNew<Input2*>(this->numtimesteps);
      xMemCpy(old_inputs,this->inputs,this->numtimesteps);
      xDelete<Input2*>(this->inputs);
   }
 
   this->numtimesteps=this->numtimesteps+1;
   this->timesteps=xNew<IssmDouble>(this->numtimesteps);
   this->inputs   = xNew<Input2*>(this->numtimesteps);
 
   if (this->numtimesteps > 1){
      xMemCpy(this->inputs,old_inputs,this->numtimesteps-1);
      xMemCpy(this->timesteps,old_timesteps,this->numtimesteps-1);
      xDelete(old_timesteps);
      xDelete<Input2*>(old_inputs);
   }
 
   /*go ahead and plug: */
   this->timesteps[this->numtimesteps-1] = time;
   this->inputs[this->numtimesteps-1]    = NULL;
   this->AddPentaTimeInput(this->numtimesteps-1,numindices,indices,values_in,interp_in);
 
}
/*}}}*/
void TransientInput2::AddTriaTimeInput(int step,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/
 
   _assert_(step>=0 && step<this->numtimesteps);
 
   /*Create it if necessary*/
   if(this->inputs[step]){
      if(this->inputs[step]->ObjectEnum()!=TriaInput2Enum) _error_("cannot add Element values to a "<<EnumToStringx(this->inputs[step]->ObjectEnum()));
   }
   else{
      this->inputs[step] = new TriaInput2(this->numberofelements_local,this->numberofvertices_local,interp_in);
   }
 
   /*Set input*/
   TriaInput2* input = xDynamicCast<TriaInput2*>(this->inputs[step]);
   input->SetInput(interp_in,numindices,indices,values_in);
 
}
/*}}}*/
void TransientInput2::AddPentaTimeInput(int step,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/
 
   _assert_(step>=0 && step<this->numtimesteps);
 
   /*Create it if necessary*/
   if(this->inputs[step]){
      if(this->inputs[step]->ObjectEnum()!=PentaInput2Enum) _error_("cannot add Element values to a "<<EnumToStringx(this->inputs[step]->ObjectEnum()));
   }
   else{
      this->inputs[step] = new PentaInput2(this->numberofelements_local,this->numberofvertices_local,interp_in);
   }
 
   /*Set input*/
   PentaInput2* input = xDynamicCast<PentaInput2*>(this->inputs[step]);
   input->SetInput(interp_in,numindices,indices,values_in);
 
}
/*}}}*/
void TransientInput2::GetAllTimes(IssmDouble** ptimesteps,int* pnumtimesteps){/*{{{*/
 
   if(ptimesteps){
      *ptimesteps=xNew<IssmDouble>(this->numtimesteps);
      xMemCpy(*ptimesteps,this->timesteps,this->numtimesteps);
   }
   if(pnumtimesteps){
      *pnumtimesteps = this->numtimesteps;
   }
 
}
/*}}}*/
TriaInput2* TransientInput2::GetTriaInput(){/*{{{*/
 
   IssmDouble time;
   this->parameters->FindParam(&time,TimeEnum);
   return this->GetTriaInput(time);
 
}
/*}}}*/
TriaInput2* TransientInput2::GetTriaInput(IssmDouble time){/*{{{*/
 
   /*Set current time input*/
   this->SetCurrentTimeInput(time);
   _assert_(this->current_input);
 
   /*Cast and return*/
   if(this->current_input->ObjectEnum()!=TriaInput2Enum){
      _error_("Cannot return a TriaInput2");
   }
   return xDynamicCast<TriaInput2*>(this->current_input);
 
}
/*}}}*/
TriaInput2* TransientInput2::GetTriaInput(IssmDouble start_time, IssmDouble end_time, int averaging_method){/*{{{*/
 
   /*Set current time input*/
   this->SetAverageAsCurrentTimeInput(start_time,end_time,averaging_method);
   _assert_(this->current_input);
 
   /*Cast and return*/
   if(this->current_input->ObjectEnum()!=TriaInput2Enum){
      _error_("Cannot return a TriaInput2");
   }
   return xDynamicCast<TriaInput2*>(this->current_input);
 
}
/*}}}*/
TriaInput2* TransientInput2::GetTriaInput(int offset){/*{{{*/
 
   /*Check offset*/
   if(offset<0 || offset>this->numtimesteps-1){
      _error_("Cannot return input for offset "<<offset);
   }
   Input2* input = this->inputs[offset];
 
   /*Cast and return*/
   _assert_(input);
   if(input->ObjectEnum()!=TriaInput2Enum) _error_("Cannot return a TriaInput2");
   return xDynamicCast<TriaInput2*>(input);
 
}
/*}}}*/
PentaInput2* TransientInput2::GetPentaInput(){/*{{{*/
 
   IssmDouble time;
   this->parameters->FindParam(&time,TimeEnum);
   return this->GetPentaInput(time);
}
/*}}}*/
PentaInput2* TransientInput2::GetPentaInput(IssmDouble time){/*{{{*/
 
   /*Set current time input*/
   this->SetCurrentTimeInput(time);
   _assert_(this->current_input);
 
   /*Cast and return*/
   if(this->current_input->ObjectEnum()!=PentaInput2Enum){
      _error_("Cannot return a PentaInput2");
   }
   return xDynamicCast<PentaInput2*>(this->current_input);
 
}
/*}}}*/
PentaInput2* TransientInput2::GetPentaInput(int offset){/*{{{*/
 
 
   /*Check offset*/
   if(offset<0 || offset>this->numtimesteps-1){
      _error_("Cannot return input for offset "<<offset);
   }
   Input2* input = this->inputs[offset];
 
   /*Cast and return*/
   if(input->ObjectEnum()!=PentaInput2Enum) _error_("Cannot return a PentaInput2");
   return xDynamicCast<PentaInput2*>(input);
 
}
/*}}}*/
PentaInput2* TransientInput2::GetPentaInput(IssmDouble start_time, IssmDouble end_time, int averaging_method){/*{{{*/
 
   /*Set current time input*/
   this->SetAverageAsCurrentTimeInput(start_time,end_time,averaging_method);
   _assert_(this->current_input);
 
   /*Cast and return*/
   if(this->current_input->ObjectEnum()!=PentaInput2Enum){
      _error_("Cannot return a PentaInput2");
   }
   return xDynamicCast<PentaInput2*>(this->current_input);
 
}
/*}}}*/
 
void TransientInput2::SetCurrentTimeInput(IssmDouble time){/*{{{*/
 
   /*First, recover current time from parameters: */
   bool linear_interp;
   this->parameters->FindParam(&linear_interp,TimesteppingInterpForcingsEnum);
 
   /*Figure step out*/
   int offset;
   if(!binary_search(&offset,time,this->timesteps,this->numtimesteps)){
      _error_("Input not found (is TransientInput sorted ?)");
   }
 
   if (offset==-1){
 
      /*get values for the first time: */
      _assert_(time<this->timesteps[0]);
 
      /*If already processed return*/
      if(this->current_step==0.) return;
 
      /*Prepare input*/
      if(this->current_input) delete this->current_input;
      this->current_step = 0.;
      this->current_input = this->inputs[0]->copy();
 
   }
   else if(offset==(this->numtimesteps-1) || !linear_interp){
 
      /*get values for the last time: */
      _assert_(time>=this->timesteps[offset]);
 
      /*If already processed return*/
      if(this->current_step==reCast<IssmDouble>(offset)) return;
 
      /*Prepare input*/
      if(this->current_input) delete this->current_input;
      this->current_step  = reCast<IssmDouble>(offset);
      this->current_input = this->inputs[offset]->copy();
   }
   else {
 
      /*Interpolate */
      _assert_(time>=this->timesteps[offset] && time<this->timesteps[offset+1]);
 
      /*get values between two times [offset:offset+1[, Interpolate linearly*/
      IssmDouble deltat=this->timesteps[offset+1]-this->timesteps[offset];
      IssmDouble this_step = reCast<IssmDouble>(offset) + (time - this->timesteps[offset])/deltat;
 
      /*If already processed return*/
      if(fabs(this->current_step-this_step)<1.e-5) return;
      //    if(this->current_step>this_step-1.e-5 && this->current_step<this_step+1.e-5) return;
 
      /*Prepare input*/
      if(this->current_input) delete this->current_input;
      this->current_step = this_step;
      IssmDouble alpha2=(time-this->timesteps[offset])/deltat;
      IssmDouble alpha1=(1.0-alpha2);
 
      Input2* input1=this->inputs[offset];
      Input2* input2=this->inputs[offset+1];
 
      this->current_input = input1->copy();
      this->current_input->Scale(alpha1);
      this->current_input->AXPY(input2,alpha2);
   }
 
}/*}}}*/
void TransientInput2::SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time, int averaging_method){/*{{{*/
 
   IssmDouble  dt,durinv;
   IssmDouble  dtsum=0;
   IssmDouble  timespan,mid_step;
   int         found,start_offset,end_offset,input_offset;
 
   /*go through the timesteps, and grab offset for start and end*/
   found=binary_search(&start_offset,start_time,this->timesteps,this->numtimesteps);
   if(!found) _error_("Input not found (is TransientInput sorted ?)");
   found=binary_search(&end_offset,end_time,this->timesteps,this->numtimesteps);
   if(!found) _error_("Input not found (is TransientInput sorted ?)");
 
   if(start_offset==-1){
      timespan=this->timesteps[end_offset]-start_time;
   }
   else{
      timespan=this->timesteps[end_offset]-this->timesteps[start_offset];
   }
   mid_step=reCast<IssmDouble>(start_offset)+0.5*timespan;
 
   /*If already processed return, we set step in the middle of the interval*/
   if(fabs(this->current_step-mid_step)<1.e-5) return;
   /*If not processed set current_step*/
   if(this->current_input) delete this->current_input;
   this->current_step = mid_step;
 
   int offset=start_offset;
   while(offset < end_offset){
      if(offset==start_offset){
         dt=this->timesteps[offset+1]-start_time;
         _assert_(dt>0.);
         if(offset==end_offset-1){
            dt=end_time-start_time;
            _assert_(dt>0.);
         }
      }
      else if(offset==end_offset-1){
         dt=end_time-this->timesteps[offset];
         _assert_(dt>0.);
      }
      else{
         dt=this->timesteps[offset+1]-this->timesteps[offset];
         _assert_(dt>0.);
      }
      Input2* stepinput=this->inputs[offset+1];
 
      switch(averaging_method){
         case 0: /*Arithmetic mean*/
            if(offset==start_offset){
               this->current_input=stepinput->copy();
               this->current_input->Scale(dt);
            }
            else{
               this->current_input->AXPY(stepinput,dt);
            }
            break;
         case 1: /*Geometric mean*/
            if(offset==start_offset){
               this->current_input = stepinput->copy();
               this->current_input->Scale(dt);
            }
            else{
               stepinput->Scale(dt);
               this->current_input->PointWiseMult(stepinput);
            }
            break;
         case 2: /*Harmonic mean*/
            if(offset==start_offset){
               this->current_input = stepinput->copy();
               this->current_input->Pow(-1);
               this->current_input->Scale(dt);
            }
            else{
               stepinput->Pow(-1);
               this->current_input->AXPY(stepinput,dt);
            }
         default:
            _error_("averaging method is not recognised");
      }
      dtsum+=dt;
      offset+=1;
   }
   _assert_(dtsum>0);
   durinv=1./dtsum;
   /*Integration done, now normalize*/
   switch(averaging_method){
      case 0: //Arithmetic mean
         this->current_input->Scale(durinv);
         break;
      case 1: /*Geometric mean*/
         this->current_input->Pow(durinv);
         break;
      case 2: /*Harmonic mean*/
         this->current_input->Scale(durinv);
         this->current_input->Pow(-1);
      default:
         _error_("averaging method is not recognised");
   }
}/*}}}*/
IssmDouble  TransientInput2::GetTimeByOffset(int offset){/*{{{*/
   if(offset<0) offset=0;
   _assert_(offset<this->numtimesteps);
   return this->timesteps[offset];
}
/*}}}*/
int  TransientInput2::GetTimeInputOffset(IssmDouble time){/*{{{*/
 
   int offset;
 
   /*go through the timesteps, and figure out which interval we
    *     *fall within. Then interpolate the values on this interval: */
   int found=binary_search(&offset,time,this->timesteps,this->numtimesteps);
   if(!found) _error_("Input not found (is TransientInput sorted ?)");
 
   return offset;
}
/*}}}*/