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

Timestamp:

08/25/22 16:50:29 (3 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 27230

Location:

issm/trunk

Files:

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

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

@@ -9 +9 @@
 #endif
 
+#include <numeric>
 #include "./TransientInput.h"
 #include "./TriaInput.h"
@@ -208 +209 @@
 
         /*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");
+        if(this->numtimesteps>0 && time<this->timesteps[this->numtimesteps-1]) _error_("timestep values must increase sequentially, here " << this->timesteps[this->numtimesteps-1] <<" is the last step but smaller than the preceding "<< time<<"\n");
 
         IssmDouble *old_timesteps = NULL;
@@ -433 +434 @@
 
         /*First, recover current time from parameters: */
-        bool linear_interp,cycle;
+        bool linear_interp,average,cycle;
         IssmDouble dt;
         this->parameters->FindParam(&linear_interp,TimesteppingInterpForcingEnum);
+        this->parameters->FindParam(&average,TimesteppingAverageForcingEnum);
         this->parameters->FindParam(&cycle,TimesteppingCycleForcingEnum);
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);          /*transient core time step*/
@@ -471 +473 @@
 
         /*Figure step out*/
-        int offset;
+        int offset, prevoffset;
         if(!binary_search(&offset,time,this->timesteps,this->numtimesteps)){
+                _error_("Input not found (is TransientInput sorted ?)");
+        }
+        if(!binary_search(&prevoffset,reCast<IssmDouble>(time-dt),this->timesteps,this->numtimesteps)){
                 _error_("Input not found (is TransientInput sorted ?)");
         }
@@ -488 +493 @@
                 this->current_step = 0.;
                 this->current_input = this->inputs[0]->copy();
+
+        }
+        else if(offset-prevoffset>1 && prevoffset >=0 && average){
+                /*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[prevoffset]->copy();
+                for(int i=prevoffset+1;i<offset;i++) {
+                        this->current_input->AXPY(this->inputs[i],+1.0);
+                }
+                this->current_input->Scale(1./(offset-prevoffset));
 
         }
@@ -535 +558 @@
         IssmDouble  timespan,mid_step;
         int         found,start_offset,end_offset,input_offset;
+
 
         /*go through the timesteps, and grab offset for start and end*/
@@ -573 +597 @@
                         _assert_(dt>0.);
                 }
-                Input* stepinput=this->inputs[offset+1];
-
+                Input* stepinput=this->inputs[offset+1]->copy();
                 switch(averaging_method){
                         case 0: /*Arithmetic mean*/
@@ -608 +631 @@
                                 _error_("averaging method is not recognised");
                 }
+                delete stepinput;
                 dtsum+=dt;
                 offset+=1;