Changeset 24313 for issm/trunk/src/c/classes/Inputs/TransientInput.cpp

Timestamp:

11/01/19 12:01:57 (5 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 24310

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/c (modified) (1 prop)
• src/c/classes/Inputs/TransientInput.cpp (modified) (9 diffs)

issm/trunk/src/c

@@ -5 +5 @@
 *.obj
 *.exe
+*.mod
 appscan.*
 ouncemake*
@@ -14 +15 @@
 probe.results
 stXXXX*
-*.deps
-*.dirstamp
+.deps
+.dirstamp
 .libs
 issm
@@ -21 +22 @@
 issm_slr
 issm_ocean
-lnb_param.mod
-lovenb_sub.mod
-model.mod
-util.mod
+issm_dakota

@@ -5 +5 @@
 *.obj
 *.exe
+*.mod
 appscan.*
 ouncemake*
@@ -14 +15 @@
 probe.results
 stXXXX*
-*.deps
-*.dirstamp
+.deps
+.dirstamp
 .libs
 issm
@@ -21 +22 @@
 issm_slr
 issm_ocean
-lnb_param.mod
-lovenb_sub.mod
-model.mod
-util.mod
+issm_dakota

issm/trunk/src/c/classes/Inputs/TransientInput.cpp

@@ -60 +60 @@
 }
 /*}}}*/
+void TransientInput::AddTimeInput(Input* input,IssmDouble time){/*{{{*/
+
+        /*insert values at time step: */
+        if (this->numtimesteps>0 && time<=this->timesteps[this->numtimesteps-1]) _error_("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<IssmDouble>(this->numtimesteps);
+                xMemCpy(old_timesteps,this->timesteps,this->numtimesteps);
+                xDelete(this->timesteps);
+        }
+
+        this->numtimesteps=this->numtimesteps+1;
+        this->timesteps=xNew<IssmDouble>(this->numtimesteps);
+
+        if (this->numtimesteps > 1){
+                xMemCpy(this->timesteps,old_timesteps,this->numtimesteps-1);
+                xDelete(old_timesteps);
+        }
+
+        /*go ahead and plug: */
+        this->timesteps[this->numtimesteps-1]=time;
+        inputs->AddObject(input);
+
+}
+/*}}}*/
+void TransientInput::AddTimeInput(Input* input){/*{{{*/
+
+        _assert_(this->inputs->Size()<this->numtimesteps);
+        inputs->AddObject(input);
+
+}
+/*}}}*/
 
 /*Object virtual functions definitions:*/
@@ -87 +122 @@
         _printf_("   enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
         _printf_("   numtimesteps: " << this->numtimesteps << "\n");
-        _printf_("---inputs: \n"); 
+        _printf_("---inputs: \n");
         for(i=0;i<this->numtimesteps;i++){
                 _printf_("   time: " << this->timesteps[i]<<"  ");
@@ -120 +155 @@
 
 /*TransientInput management*/
+void TransientInput::Configure(Parameters* parameters){/*{{{*/
+        this->parameters=parameters;
+}
+/*}}}*/
+int  TransientInput::GetResultArraySize(void){/*{{{*/
+
+        return 1;
+}
+/*}}}*/
+int  TransientInput::GetResultInterpolation(void){/*{{{*/
+
+        IssmDouble time;
+        int        output;
+
+        parameters->FindParam(&time,TimeEnum);
+        Input* input=GetTimeInput(time);
+        output = input->GetResultInterpolation();
+
+        /*Clean up and return*/
+        delete input;
+        return output;
+
+}
+/*}}}*/
+int  TransientInput::GetResultNumberOfNodes(void){/*{{{*/
+
+        IssmDouble time;
+        int        output;
+
+        parameters->FindParam(&time,TimeEnum);
+        Input* input=GetTimeInput(time);
+        output = input->GetResultNumberOfNodes();
+
+        /*Clean up and return*/
+        delete input;
+        return output;
+
+}
+/*}}}*/
 int TransientInput::InstanceEnum(void){/*{{{*/
 
@@ -174 +248 @@
 }
 /*}}}*/
-void TransientInput::GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){/*{{{*/
+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;
+
+}
+/*}}}*/
+void TransientInput::GetInputAverageOverTimeSlice(IssmDouble* pvalue, Gauss* gauss, IssmDouble start_time,IssmDouble end_time){/*{{{*/
+
+        int averaging = 0;
+
+        Input* input=GetInputAverageOverTime(start_time,end_time,averaging);
+
+        /*Call input function*/
+        input->GetInputValue(pvalue, gauss);
+
+        //*pvalues=values;
+        delete input;
+}
+/*}}}*/
+Input* TransientInput::GetInputAverageOverTime(IssmDouble start_time,IssmDouble end_time,int averaging){/*{{{*/
+
+        int found;
+        int offset,start_offset,end_offset;
+        IssmDouble subdt;
+        IssmDouble slice_duration;
+
+        this->parameters->FindParam(&subdt,TimesteppingTimeStepEnum); //duration of each substeps
+
+        Input *averageinput  = NULL;
+
+
+        slice_duration=end_time-start_time;
+        start_time+=subdt; //because time is actually considered at the end of the timestep
+
+        /*go through the timesteps, and grab offset for the start of the slice*/
+        found=binary_search(&offset,start_time,this->timesteps,this->numtimesteps);
+        if(!found) _error_("Input not found (is TransientInput sorted ?)");
+        /*go through the timesteps, and grab offset for the end of the slice*/
+        found=binary_search(&end_offset,end_time,this->timesteps,this->numtimesteps);
+        if(!found) _error_("Input not found (is TransientInput sorted ?)");
+
+        start_offset = offset;
+        //stack the input for each timestep in the slice
+        while(offset <= end_offset ){
+                Input *currentinput  = NULL;
+
+                if (offset==-1){
+                        /*get values for the first time: */
+                        _assert_(start_time<this->timesteps[0]);
+                        currentinput=(Input*)((Input*)this->inputs->GetObjectByOffset(0))->copy();
+                }
+                else if(offset==(this->numtimesteps-1)){
+                        /*get values for the last time: */
+                        _assert_(end_time>=this->timesteps[offset]);
+                        currentinput=(Input*)((Input*)this->inputs->GetObjectByOffset(offset))->copy();
+                }
+                else{
+                        currentinput=(Input*)((Input*)this->inputs->GetObjectByOffset(offset))->copy();
+                }
+                switch(averaging){
+                        case 0: //Arithmetic mean
+                                if (offset==start_offset){
+                                        averageinput=(Input*)currentinput->copy();
+                                        averageinput->Scale(subdt);
+                                }
+                                else{
+                                        averageinput->AXPY(currentinput,subdt);
+                                }
+                                break;
+                        case 1: //Geometric mean
+                                if (offset==start_offset){
+                                        averageinput=(Input*)currentinput->copy();
+                                        averageinput->Scale(subdt);
+                                }
+                                else{
+                                        currentinput->Scale(subdt);
+                                        averageinput->PointwiseMult(currentinput);
+                                }
+                                break;
+                        case 2: //Harmonic mean
+                                if (offset==start_offset){
+                                        averageinput=(Input*)currentinput->copy();
+                                        averageinput->Pow(-1);
+                                        averageinput->Scale(subdt);
+                                }
+                                else{
+                                        currentinput->Pow(-1);
+                                        averageinput->AXPY(currentinput,subdt);
+                                }
+                                break;
+                        default:
+                                _error_("averaging method is not recognised");
+                }
+                delete currentinput;
+                offset+=1;
+        }
+
+        //summation is done, now we normalise
+        switch(averaging){
+                case 0: //Arithmetic mean
+                        averageinput->Scale(1.0/slice_duration);
+                        break;
+                case 1: //Geometric mean
+                        averageinput->Pow(1.0/slice_duration);
+                        break;
+                case 2: //Harmonic mean
+                        averageinput->Scale(1.0/slice_duration);
+                        averageinput->Pow(-1.0);
+                        break;
+                default:
+                        _error_("averaging method is not recognised");
+        }
+        return averageinput;
+}
+/*}}}*/
+void TransientInput::GetInputAveragesOnAllTime(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){/*{{{*/
 
         int i;
@@ -195 +395 @@
 }
 /*}}}*/
-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;
-
-}
-/*}}}*/
-void TransientInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, Gauss* 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;
-}
-/*}}}*/
-void TransientInput::GetInputValue(IssmDouble* pvalue,Gauss* 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->GetInputValue(pvalue,gauss);
-
-        delete input;
-}
-/*}}}*/
-void TransientInput::GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){/*{{{*/
-
-        /*Retrieve interpolated values for this time step: */
-        Input* input=GetTimeInput(time);
-
-        /*Call input function*/
-        input->GetInputValue(pvalue,gauss);
-
-        delete input;
-}
-/*}}}*/
-int  TransientInput::GetTimeInputOffset(IssmDouble time){/*{{{*/
-
-        int        found;
-        int        offset;
-
-        /*go through the timesteps, and figure out which interval we
-         *     *fall within. Then interpolate the values on this interval: */
-        found=binary_search(&offset,time,this->timesteps,this->numtimesteps);
-        if(!found) _error_("Input not found (is TransientInput sorted ?)");
-
-        return offset;
-}
-/*}}}*/
-IssmDouble  TransientInput::GetTimeByOffset(int offset){/*{{{*/
-        if (offset < 0) offset=0;
-        _assert_(offset<(this->numtimesteps));
-        return this->timesteps[offset];
-}
-/*}}}*/
-void TransientInput::GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){/*{{{*/
+void TransientInput::GetInputAveragesUpToCurrentTime(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){/*{{{*/
 
         int         i;
@@ -313 +435 @@
 }
 /*}}}*/
-void TransientInput::GetInputAverageOnTimes(IssmDouble** pvalues, IssmDouble init_time){/*{{{*/
-
-        int                                     numnodes;
-        IssmDouble  subtimestep;
-        IssmDouble* values= NULL;
-        Input* input=NULL;
-
-        /*Initialize numnode from transientinput out of time loop*/
-        for(int i=0;i<this->numtimesteps;i++){
-                /*First compute the lengt of the sub-timestep*/
-                if(i==0){
-                        subtimestep = this->timesteps[i]-init_time;
-                }
-                else{
-                        subtimestep = this->timesteps[i]-this->timesteps[i-1];
-                }       
-                /*Figure out type of input and get it*/
-                input = xDynamicCast<Input*>(this->inputs->GetObjectByOffset(i)); _assert_(input);
-                switch(input->ObjectEnum()){
-                case TriaInputEnum:{
-                        TriaInput* triainput = (TriaInput*)this->inputs->GetObjectByOffset(i); _assert_(triainput);
-                        if(i==0){
-                                numnodes= triainput->NumberofNodes(triainput->interpolation_type);
-                                values=xNewZeroInit<IssmDouble>(numnodes);
-                        }
-                        /*Sum the values weighted by their respective sub-timestep*/
-                        for(int j=0;j<numnodes;j++){
-                                values[j]+=(triainput->values[j]*subtimestep);
-                        }
-                        break;
-                }
-                case PentaInputEnum:{
-                        PentaInput*     pentainput = (PentaInput*)this->inputs->GetObjectByOffset(i); _assert_(pentainput);
-                        if(i==0){
-                                numnodes= pentainput->NumberofNodes(pentainput->interpolation_type);
-                                values=xNewZeroInit<IssmDouble>(numnodes);
-                        }
-                        /*Sum the values weighted by their respective sub-timestep*/
-                        for(int j=0;j<numnodes;j++){
-                                values[j]+=(pentainput->values[j]*subtimestep);
-                        }
-                        break;
-                }
-                case TetraInputEnum:{
-                        TetraInput*     tetrainput = (TetraInput*)this->inputs->GetObjectByOffset(i); _assert_(tetrainput);
-                        if(i==0){
-                                numnodes= tetrainput->NumberofNodes(tetrainput->interpolation_type);
-                                values=xNewZeroInit<IssmDouble>(numnodes);
-                        }
-                        /*Sum the values weighted by their respective sub-timestep*/
-                        for(int j=0;j<numnodes;j++){
-                                values[j]+=(tetrainput->values[j]*subtimestep);
-                        }
-                        break;
-                }
-                default: _error_("not implemented yet");
-                }
-        }
-        /*Compute the average based on the length of the full timestep*/
-        for(int j=0;j<numnodes;j++){
-                values[j]=values[j]/(this->timesteps[this->numtimesteps-1]-init_time);
-        }
-        *pvalues=values;
-}
-/*}}}*/
-
+void TransientInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, Gauss* 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;
+}
+/*}}}*/
+void TransientInput::GetInputValue(IssmDouble* pvalue,Gauss* 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->GetInputValue(pvalue,gauss);
+
+        delete input;
+}
+/*}}}*/
+void TransientInput::GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){/*{{{*/
+
+        /*Retrieve interpolated values for this time step: */
+        Input* input=GetTimeInput(time);
+
+        /*Call input function*/
+        input->GetInputValue(pvalue,gauss);
+
+        delete input;
+}
+/*}}}*/
+IssmDouble  TransientInput::GetTimeByOffset(int offset){/*{{{*/
+        if (offset < 0) offset=0;
+        _assert_(offset<(this->numtimesteps));
+        return this->timesteps[offset];
+}
+/*}}}*/
+int  TransientInput::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;
+}
+/*}}}*/
 /*Intermediary*/
-void TransientInput::AddTimeInput(Input* input,IssmDouble time){/*{{{*/
-
-        /*insert values at time step: */
-        if (this->numtimesteps>0 && time<=this->timesteps[this->numtimesteps-1]) _error_("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<IssmDouble>(this->numtimesteps);
-                xMemCpy(old_timesteps,this->timesteps,this->numtimesteps);
-                xDelete(this->timesteps);
-        }
-
-        this->numtimesteps=this->numtimesteps+1;
-        this->timesteps=xNew<IssmDouble>(this->numtimesteps);
-
-        if (this->numtimesteps > 1){
-                xMemCpy(this->timesteps,old_timesteps,this->numtimesteps-1);
-                xDelete(old_timesteps);
-        }
-
-        /*go ahead and plug: */
-        this->timesteps[this->numtimesteps-1]=time;
-        inputs->AddObject(input);
-
-}
-/*}}}*/
-void TransientInput::AddTimeInput(Input* input){/*{{{*/
-
-        _assert_(this->inputs->Size()<this->numtimesteps);
-        inputs->AddObject(input);
-
-}
-/*}}}*/
-void TransientInput::Configure(Parameters* parameters){/*{{{*/
-        this->parameters=parameters;
-}
-/*}}}*/
 void TransientInput::Extrude(int start){/*{{{*/
 
@@ -424 +501 @@
                 ((Input*)this->inputs->GetObjectByOffset(i))->Extrude(start);
         }
-}
-/*}}}*/
-int  TransientInput::GetResultArraySize(void){/*{{{*/
-
-        return 1;
-}
-/*}}}*/
-int  TransientInput::GetResultInterpolation(void){/*{{{*/
-
-        IssmDouble time;
-        int        output;
-
-        parameters->FindParam(&time,TimeEnum);
-        Input* input=GetTimeInput(time);
-        output = input->GetResultInterpolation();
-
-        /*Clean up and return*/
-        delete input;
-        return output;
-
-}
-/*}}}*/
-int  TransientInput::GetResultNumberOfNodes(void){/*{{{*/
-
-        IssmDouble time;
-        int        output;
-
-        parameters->FindParam(&time,TimeEnum);
-        Input* input=GetTimeInput(time);
-        output = input->GetResultNumberOfNodes();
-
-        /*Clean up and return*/
-        delete input;
-        return output;
-
 }
 /*}}}*/
@@ -476 +518 @@
         Input *input2 = NULL;
 
-        /*go through the timesteps, and figure out which interval we 
+        /*go through the timesteps, and figure out which interval we
          *fall within. Then interpolate the values on this interval: */
         found=binary_search(&offset,intime,this->timesteps,this->numtimesteps);
@@ -498 +540 @@
                 alpha1=(1.0-alpha2);
 
-                input1=(Input*)this->inputs->GetObjectByOffset(offset); 
+                input1=(Input*)this->inputs->GetObjectByOffset(offset);
                 input2=(Input*)this->inputs->GetObjectByOffset(offset+1);
 

@@ -5 +5 @@
 *.obj
 *.exe
+*.mod
 appscan.*
 ouncemake*
@@ -14 +15 @@
 probe.results
 stXXXX*
-*.deps
-*.dirstamp
+.deps
+.dirstamp
 .libs
 issm
@@ -21 +22 @@
 issm_slr
 issm_ocean
-lnb_param.mod
-lovenb_sub.mod
-model.mod
-util.mod
+issm_dakota