Changeset 25024

Timestamp:

06/12/20 08:22:16 (5 years ago)

Author:

bdef

Message:

CHG:moving averaging on inputs and modification on 905 test

Location:

issm/trunk-jpl

Files:

• src/c/analyses/HydrologyDCEfficientAnalysis.cpp (modified) (4 diffs)
• src/c/analyses/HydrologyDCInefficientAnalysis.cpp (modified) (2 diffs)
• src/c/classes/Elements/Penta.cpp (modified) (1 diff)
• src/c/classes/Elements/Tria.cpp (modified) (24 diffs)
• src/c/classes/Inputs2/Inputs2.cpp (modified) (1 diff)
• src/c/classes/Inputs2/Inputs2.h (modified) (1 diff)
• src/c/classes/Inputs2/TransientInput2.cpp (modified) (4 diffs)
• src/c/classes/Inputs2/TransientInput2.h (modified) (1 diff)
• src/c/cores/hydrology_core.cpp (modified) (2 diffs)
• src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp (modified) (1 diff)
• src/m/classes/hydrologydc.m (modified) (1 diff)
• test/Archives/Archive905.arch (modified) ( previous)
• test/NightlyRun/test905.m (modified) (2 diffs)
• test/NightlyRun/test905.py (modified) (4 diffs)

issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

@@ -299 +299 @@
         /*Intermediaries */
         int        smb_model;
+        int        smb_averaging;
+        int        smbsubstepping;
+        int        hydrologysubstepping;
         IssmDouble dt,scalar,water_head;
         IssmDouble water_load,transfer,runoff_value;
@@ -307 +310 @@
         IssmDouble *xyz_list             = NULL;
         Input2     *old_wh_input         = NULL;
+        Input2     *dummy_input          = NULL;
         Input2     *surface_runoff_input = NULL;
 
@@ -321 +325 @@
         basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
         basalelement ->FindParam(&smb_model,SmbEnum);
+        basalelement->FindParam(&smb_averaging,SmbAveragingEnum);
 
         Input2* epl_thick_input  = basalelement->GetInput2(HydrologydcEplThicknessSubstepEnum); _assert_(epl_thick_input);
@@ -329 +334 @@
         Input2* base_input       = basalelement->GetInput2(BaseEnum); _assert_(base_input);
 
+        IssmDouble time;
+        basalelement->FindParam(&time,TimeEnum);
+
         if(dt!= 0.){
                 old_wh_input = basalelement->GetInput2(EplHeadOldEnum);            _assert_(old_wh_input);
         }
         if(smb_model==SMBgradientscomponentsEnum){
-                surface_runoff_input = basalelement->GetInput2(SmbRunoffEnum); _assert_(surface_runoff_input);
+                basalelement->FindParam(&smbsubstepping,SmbStepsPerStepEnum);
+                basalelement->FindParam(&hydrologysubstepping,HydrologyStepsPerStepEnum);
+
+                if(smbsubstepping==1){
+                        //no substeping for the smb we take the result from there
+                        dummy_input = basalelement->GetInput2(SmbRunoffEnum); _assert_(dummy_input);
+                }
+                else if(smbsubstepping>1 && smbsubstepping<=hydrologysubstepping){
+                        //finer hydro stepping, we take the value at the needed time
+                        dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum, time); _assert_(dummy_input);
+                }
+                else{
+                        //finer stepping in smb, we average the runoff from transient input
+                        dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time,smb_averaging); _assert_(dummy_input);
+                }
+                surface_runoff_input=xDynamicCast<Input2*>(dummy_input); _assert_(surface_runoff_input);
         }
 

issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

@@ -375 +375 @@
         basalelement->FindParam(&time,TimeEnum);
 
-
         if(dt!= 0.){
                 old_wh_input = basalelement->GetInput2(SedimentHeadOldEnum); _assert_(old_wh_input);
@@ -384 +383 @@
 
                 if(smbsubstepping==1){
+                        //no substeping for the smb we take the result from there
                         dummy_input = basalelement->GetInput2(SmbRunoffEnum); _assert_(dummy_input);
                 }
                 else if(smbsubstepping>1 && smbsubstepping<=hydrologysubstepping){
+                        //finer hydro stepping, we take the value at the needed time
                         dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum, time); _assert_(dummy_input);
                 }
                 else{
+                        //finer stepping in smb, we average the runoff from transient input
                         dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time,smb_averaging); _assert_(dummy_input);
                 }

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -1009 +1009 @@
 /*}}}*/
 void       Penta::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time,int averaging_method){/*{{{*/
-
         _assert_(end_time>start_time);
 
-        /*Intermediaries*/
-        IssmDouble averaged_values[NUMVERTICES];
-        IssmDouble current_values[NUMVERTICES];
-        IssmDouble dt;
-        int        found,start_offset,end_offset;
-
         /*Get transient input time steps*/
-        int         numtimesteps;
-        IssmDouble *timesteps    = NULL;
         TransientInput2* transient_input  = this->inputs2->GetTransientInput(transientinput_enum);
-        transient_input->GetAllTimes(&timesteps,&numtimesteps);
-
-        /*go through the timesteps, and grab offset for start and end*/
-        found=binary_search(&start_offset,start_time,timesteps,numtimesteps);
-        if(!found) _error_("Input not found (is TransientInput sorted ?)");
-        found=binary_search(&end_offset,end_time,timesteps,numtimesteps);
-        if(!found) _error_("Input not found (is TransientInput sorted ?)");
-
-        Gauss* gauss=this->NewGauss();
-
-        /*stack the input for each timestep in the slice*/
-        int offset = start_offset;
-        while(offset < end_offset ){
-
-
-                if(offset==-1){
-                        /*get values for the first time: */
-                        _assert_(start_time<timesteps[0]);
-                        PentaInput2* input = transient_input->GetPentaInput(0);
-
-                        int interpolation = input->GetInterpolation();
-                        _assert_(interpolation==P1Enum);
-                        /*Intermediaries*/
-                        int numindices;
-                        int indices[NUMVERTICES];
-                        numindices = NUMVERTICES;
-                        for(int i=0;i<NUMVERTICES;i++) indices[i] = vertices[i]->lid;
-                        input->Serve(numindices,&indices[0]);
-
-                        for(int iv=0;iv<NUMVERTICES;iv++){
-                                gauss->GaussVertex(iv);
-                                input->GetInputValue(&current_values[iv],gauss);
-                        }
-                        dt = timesteps[0] - start_time; _assert_(dt>0.);
-                }
-                else{
-                        PentaInput2* input = transient_input->GetPentaInput(offset+1);
-                        int interpolation = input->GetInterpolation();
-                        _assert_(interpolation==P1Enum);
-                        /*Intermediaries*/
-                        int numindices;
-                        int indices[NUMVERTICES];
-                        numindices = NUMVERTICES;
-                        for(int i=0;i<NUMVERTICES;i++) indices[i] = vertices[i]->lid;
-                        input->Serve(numindices,&indices[0]);
-
-                        for(int iv=0;iv<NUMVERTICES;iv++){
-                                gauss->GaussVertex(iv);
-                                input->GetInputValue(&current_values[iv],gauss);
-                        }
-                        if(offset == numtimesteps-1){
-                                dt = end_time - timesteps[offset]; _assert_(dt>0.);
-                        }
-                        else{
-                                dt = timesteps[offset+1] - timesteps[offset]; _assert_(dt>0.);
-                        }
-                }
-
-                switch(averaging_method){
-                        case 0: /*Arithmetic mean*/
-                                if(offset==start_offset){
-                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv]  = dt*current_values[iv];
-                                }
-                                else{
-                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] += dt*current_values[iv];
-                                }
-                                break;
-                        case 1: /*Geometric mean*/
-                                if(offset==start_offset){
-                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv]  = dt*current_values[iv];
-                                }
-                                else{
-                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] *= dt*current_values[iv];
-                                }
-                                break;
-                        case 2: /*Harmonic mean*/
-                                if(offset==start_offset){
-                                        for(int iv=0;iv<NUMVERTICES;iv++){
-                                                _assert_(current_values[iv]>1.e-50);
-                                                averaged_values[iv]  = dt*1./current_values[iv];
-                                        }
-                                }
-                                else{
-                                        for(int iv=0;iv<NUMVERTICES;iv++){
-                                                _assert_(current_values[iv]>1.e-50);
-                                                averaged_values[iv]  += dt*1./current_values[iv];
-                                        }
-                                }
-                                break;
-                        default:
-                                _error_("averaging method is not recognised");
-                }
-
-                offset+=1;
-        }
-
-        /*Integration done, now normalize*/
-        switch(averaging_method){
-                case 0: //Arithmetic mean
-                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] =  averaged_values[iv]/(end_time - start_time);
-                        break;
-                case 1: /*Geometric mean*/
-                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = pow(averaged_values[iv], 1./(end_time - start_time));
-                        break;
-                case 2: /*Harmonic mean*/
-                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = 1./(averaged_values[iv]/(end_time - start_time));
-                        break;
-                default:
-                        _error_("averaging method is not recognised");
-        }
-
-        this->AddInput2(averagedinput_enum,&averaged_values[0],P1Enum);
-
-        /*Cleanup*/
-        delete gauss;
-        xDelete<IssmDouble>(timesteps);
+        PentaInput2* averaged_input = transient_input->GetPentaInput(start_time,end_time,averaging_method);
+        Input2* averaged_copy = averaged_input->copy();
+
+        averaged_input->ChangeEnum(averagedinput_enum);
+        this->inputs2->AddInput(averaged_copy);
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -914 +914 @@
 }
 /*}}}*/
-void                        Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/
+void         Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/
 
         IssmDouble  xyz_list[NUMVERTICES][3];
@@ -1918 +1918 @@
         /*Get Input from dataset*/
         if(this->iscollapsed){
-                _error_("Get Average input not implemented in Penta yet");
+                PentaInput2* input = this->inputs2->GetPentaInput(inputenum,start_time,end_time,averaging_method);
+                if(!input) return input;
+
+                this->InputServe(input);
+                return input;
         }
         else{
@@ -2114 +2118 @@
         _assert_(end_time>start_time);
 
-        /*Intermediaries*/
-        IssmDouble averaged_values[NUMVERTICES];
-        IssmDouble current_values[NUMVERTICES];
-        IssmDouble dt;
-        int        found,start_offset,end_offset;
-
         /*Get transient input time steps*/
-        int         numtimesteps;
-        IssmDouble *timesteps    = NULL;
         TransientInput2* transient_input  = this->inputs2->GetTransientInput(transientinput_enum);
-        transient_input->GetAllTimes(&timesteps,&numtimesteps);
-
-        /*go through the timesteps, and grab offset for start and end*/
-        found=binary_search(&start_offset,start_time,timesteps,numtimesteps);
-        if(!found) _error_("Input not found (is TransientInput sorted ?)");
-        found=binary_search(&end_offset,end_time,timesteps,numtimesteps);
-        if(!found) _error_("Input not found (is TransientInput sorted ?)");
-
-        Gauss* gauss=this->NewGauss();
-
-        /*stack the input for each timestep in the slice*/
-        int offset = start_offset;
-        while(offset < end_offset ){
-                if(offset==-1){
-                        /*get values for the first time: */
-                        _assert_(start_time<timesteps[0]);
-                        TriaInput2* input = transient_input->GetTriaInput(0);
-                        _assert_(input->GetInterpolation()==P1Enum);
-                        this->InputServe(input);
-                        for(int iv=0;iv<NUMVERTICES;iv++){
-                                gauss->GaussVertex(iv);
-                                input->GetInputValue(&current_values[iv],gauss);
-                        }
-                        dt = timesteps[0] - start_time; _assert_(dt>0.);
-                }
-                else{
-                        TriaInput2* input = transient_input->GetTriaInput(offset+1);
-                        _assert_(input->GetInterpolation()==P1Enum);
-                        this->InputServe(input);
-                        for(int iv=0;iv<NUMVERTICES;iv++){
-                                gauss->GaussVertex(iv);
-                                input->GetInputValue(&current_values[iv],gauss);
-                        }
-                        if(offset == numtimesteps-1){
-                                dt = end_time - timesteps[offset]; _assert_(dt>0.);
-                        }
-                        else{
-                                dt = timesteps[offset+1] - timesteps[offset]; _assert_(dt>0.);
-                        }
-                }
-
-                switch(averaging_method){
-                        case 0: /*Arithmetic mean*/
-                                if(offset==start_offset){
-                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv]  = dt*current_values[iv];
-                                }
-                                else{
-                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] += dt*current_values[iv];
-                                }
-                                break;
-                        case 1: /*Geometric mean*/
-                                if(offset==start_offset){
-                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv]  = dt*current_values[iv];
-                                }
-                                else{
-                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] *= dt*current_values[iv];
-                                }
-                                break;
-                        case 2: /*Harmonic mean*/
-                                if(offset==start_offset){
-                                        for(int iv=0;iv<NUMVERTICES;iv++){
-                                                _assert_(current_values[iv]>1.e-50);
-                                                averaged_values[iv]  = dt*1./current_values[iv];
-                                        }
-                                }
-                                else{
-                                        for(int iv=0;iv<NUMVERTICES;iv++){
-                                                _assert_(current_values[iv]>1.e-50);
-                                                averaged_values[iv]  += dt*1./current_values[iv];
-                                        }
-                                }
-                                break;
-                        default:
-                                _error_("averaging method is not recognised");
-                }
-
-                offset+=1;
-        }
-
-        /*Integration done, now normalize*/
-        switch(averaging_method){
-                case 0: //Arithmetic mean
-                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] =  averaged_values[iv]/(end_time - start_time);
-                        break;
-                case 1: /*Geometric mean*/
-                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = pow(averaged_values[iv], 1./(end_time - start_time));
-                        break;
-                case 2: /*Harmonic mean*/
-                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = (end_time - start_time)/averaged_values[iv];
-                        break;
-                default:
-                        _error_("averaging method is not recognised");
-        }
-
-        this->AddInput2(averagedinput_enum,&averaged_values[0],P1Enum);
-
-        /*Cleanup*/
-        delete gauss;
-        xDelete<IssmDouble>(timesteps);
+        TriaInput2* averaged_input = transient_input->GetTriaInput(start_time,end_time,averaging_method);
+        Input2* averaged_copy = averaged_input->copy();
+
+        averaged_copy->ChangeEnum(averagedinput_enum);
+        this->inputs2->AddInput(averaged_copy);
 }
 /*}}}*/
@@ -5329 +5231 @@
         rho_ice=FindParam(MaterialsRhoIceEnum);
         rho_earth=FindParam(MaterialsEarthDensityEnum);
-        
+
         /*recover earth area: */
         this->parameters->FindParam(&planetarea,SealevelPlanetAreaEnum);
@@ -5356 +5258 @@
         }
 
-        // correction at the north pole: given longitude of the North pole a definition 
+        // correction at the north pole: given longitude of the North pole a definition
         // closer to the other two vertices.
         if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
@@ -5362 +5264 @@
         if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
 
-        // correction at the north pole: given longitude of the North pole a definition 
+        // correction at the north pole: given longitude of the North pole a definition
         // closer to the other two vertices.
         if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
@@ -5569 +5471 @@
         masks->isiceonly[this->lid]=this->IsIceOnlyInElement();
         masks->isoceanin[this->lid]=this->IsOceanInElement();
-        
+
         /*are we fully floating:*/
         Input2* gr_input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input);
@@ -5595 +5497 @@
         IssmDouble x_element,y_element,z_element,x,y,z,dx,dy,dz,N_azim,E_azim;
 
-        #ifdef _HAVE_RESTRICT_ 
+        #ifdef _HAVE_RESTRICT_
         IssmDouble* __restrict__ G=NULL;
         IssmDouble* __restrict__ GU=NULL;
@@ -5616 +5518 @@
         IssmDouble* indices=NULL;
         #endif
- 
+
         /*elastic green function:*/
         int index;
@@ -5656 +5558 @@
                 GE=xNewZeroInit<IssmDouble>(gsize);
         }
-        
+
         /*compute area:*/
         area=GetAreaSpherical();
@@ -5676 +5578 @@
 
         constant=3/rho_earth*area/planetarea;
-        
+
         for(int i=0;i<gsize;i++){
                 IssmDouble alpha;
@@ -5687 +5589 @@
                 indices[i]=alpha/PI*reCast<IssmDouble,int>(M-1);
                 index=reCast<int,IssmDouble>(indices[i]);
-                
+
                 /*Rigid earth gravitational perturbation: */
                 if(computerigid){
@@ -5695 +5597 @@
                         G[i]+=G_elastic_precomputed[index];
                 }
-                G[i]=G[i]*constant; 
+                G[i]=G[i]*constant;
 
                 /*Elastic components:*/
@@ -5728 +5630 @@
 
         /*Free allocations:*/
-        #ifdef _HAVE_RESTRICT_ 
-        delete indices; 
-        delete G; 
-        delete GU; 
+        #ifdef _HAVE_RESTRICT_
+        delete indices;
+        delete G;
+        delete GU;
         if(horiz){
-                delete GN; 
+                delete GN;
                 delete GE;
         }
@@ -5752 +5654 @@
 
         int bp_compute_fingerprints= 0;
-                
+
         /*Compute bottom pressure contribution from ocean if requested:*/
         this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
@@ -5795 +5697 @@
         /*early return if we are fully floating:*/
         if (masks->isfullyfloating[this->lid]){
-                constant=0; 
+                constant=0;
                 #ifdef _ISSM_DEBUG_
                 this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
-                #endif 
+                #endif
                 *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
                 return;
@@ -5805 +5707 @@
         /*If we are here, we are on ice that is fully grounded or half-way to floating:*/
         if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on.
-    
+
         /*Inform mask: */
-        constant=1; 
+        constant=1;
         #ifdef _ISSM_DEBUG_
         this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
@@ -5820 +5722 @@
 
         /*retrieve precomputed G:*/
-        this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); 
+        this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
 
         /*Get area of element: precomputed in the sealevelrise_core_geometry:*/
@@ -5826 +5728 @@
 
         /*Compute fraction of the element that is grounded: */
-        if(notfullygrounded){ 
+        if(notfullygrounded){
                 IssmDouble xyz_list[NUMVERTICES][3];
                 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
                 phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias
-        } 
+        }
         else phi=1.0;
 
@@ -5892 +5794 @@
         /*elastic green function:*/
         IssmDouble* G=NULL;
-        
+
         /*water properties: */
         IssmDouble rho_water;
@@ -5908 +5810 @@
 
         /*retrieve precomputed G:*/
-        this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); 
+        this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
 
         /*Get area of element: precomputed in the sealevelrise_core_geometry:*/
@@ -5943 +5845 @@
         /*early return if we are not on the ocean:*/
         if (!masks->isoceanin[this->lid]){
-                constant=0; 
-                #ifdef _ISSM_DEBUG_ 
+                constant=0;
+                #ifdef _ISSM_DEBUG_
                 this->AddInput2(SealevelEustaticOceanMaskEnum,&constant,P0Enum);
                 #endif
                 return;
         }
-        constant=1; 
-        #ifdef _ISSM_DEBUG_ 
+        constant=1;
+        #ifdef _ISSM_DEBUG_
         this->AddInput2(SealevelEustaticOceanMaskEnum,&constant,P0Enum);
         #endif
@@ -5956 +5858 @@
         /*how many dofs are we working with here? */
         this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
-        
+
         /*retrieve precomputed G:*/
         this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&dummy); _assert_(dummy==gsize);
@@ -6004 +5906 @@
 
         /*recover elastic Green's functions for displacement:*/
-        this->inputs2->GetArrayPtr(SealevelriseGUEnum,this->lid,&GU,&gsize); 
+        this->inputs2->GetArrayPtr(SealevelriseGUEnum,this->lid,&GU,&gsize);
         if(horiz){
-                this->inputs2->GetArrayPtr(SealevelriseGEEnum,this->lid,&GE,&gsize); 
-                this->inputs2->GetArrayPtr(SealevelriseGNEnum,this->lid,&GN,&gsize); 
+                this->inputs2->GetArrayPtr(SealevelriseGEEnum,this->lid,&GE,&gsize);
+                this->inputs2->GetArrayPtr(SealevelriseGNEnum,this->lid,&GN,&gsize);
         }
 

issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp

@@ -388 +388 @@
         }
 }/*}}}*/
+PentaInput2* Inputs2::GetPentaInput(int enum_in,IssmDouble start_time,IssmDouble end_time,int averaging_method){/*{{{*/
+
+        /*Get input id*/
+        int id = EnumToIndex(enum_in);
+
+        /*Check that it has the right format*/
+        Input2* input = this->inputs[id];
+        if(!input) return NULL;
+
+        if(input->ObjectEnum()==TransientInput2Enum){
+                return xDynamicCast<TransientInput2*>(input)->GetPentaInput(start_time,end_time,averaging_method);
+        }
+        else{
+                _error_("Input "<<EnumToStringx(enum_in)<<" is not an TransientInput2");
+        }
+}/*}}}*/
 TransientInput2* Inputs2::GetTransientInput(int enum_in){/*{{{*/
 

issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.h

@@ -58 +58 @@
                 PentaInput2*     GetPentaInput(int enum_type);
                 PentaInput2*     GetPentaInput(int enum_type,IssmDouble time);
+                PentaInput2*     GetPentaInput(int enum_in,IssmDouble start_time,IssmDouble end_time,int averaging_method);
                 TransientInput2* GetTransientInput(int enum_type);
                 ElementInput2*   GetControlInput2Data(int enum_type,const char* data);

issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp

@@ -355 +355 @@
 }
 /*}}}*/
+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){/*{{{*/
 
@@ -404 +419 @@
 
                 /*If already processed return*/
-                if(this->current_step>this_step-1.e-5 && this->current_step<this_step+1.e-5) return;
+                if(abs(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*/
@@ -424 +440 @@
 
         IssmDouble  dt,durinv;
-        IssmPDouble eps=1.0e-6;
         IssmDouble  dtsum=0;
-        int         found,start_offset,end_offset;
+        IssmDouble  timespan,mid_step;
+        int         found,start_offset,end_offset,input_offset;
 
         /*go through the timesteps, and grab offset for start and end*/
-        IssmDouble temp = start_time-eps;
-        found=binary_search(&start_offset,temp,this->timesteps,this->numtimesteps);
+        found=binary_search(&start_offset,start_time,this->timesteps,this->numtimesteps);
         if(!found) _error_("Input not found (is TransientInput sorted ?)");
-        temp = end_time+eps;
-        found=binary_search(&end_offset,temp,this->timesteps,this->numtimesteps);
+        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(abs(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;
-        if(this->current_input) delete this->current_input;
         while(offset < end_offset){
-                if (offset==-1){
-                        dt=this->timesteps[0]-start_time;
-                        _assert_(start_time<this->timesteps[0]);
-                }
-                else if(offset==this->numtimesteps-1)dt=end_time-this->timesteps[offset];
-                else if(offset==start_offset && this->timesteps[offset]<start_time) dt=this->timesteps[offset+1]-start_time;
-                else if(offset==end_offset && this->timesteps[offset]>end_time) dt=end_time-this->timesteps[offset];
-                else dt=this->timesteps[offset+1]-this->timesteps[offset];
-                _assert_(dt>0.);
-
+                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];
 
@@ -454 +487 @@
                         case 0: /*Arithmetic mean*/
                                 if(offset==start_offset){
-                                        this->current_input = stepinput->copy();
+                                        this->current_input=stepinput->copy();
                                         this->current_input->Scale(dt);
                                 }

issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h

@@ -55 +55 @@
                 PentaInput2* GetPentaInput(IssmDouble time);
                 PentaInput2* GetPentaInput(int offset);
+                PentaInput2* GetPentaInput(IssmDouble start_time,IssmDouble end_time,int averaging_method);
                 Input2*      GetTimeInput(IssmDouble time){_error_("This should not happen!");};
                 IssmDouble   GetTimeByOffset(int offset);

issm/trunk-jpl/src/c/cores/hydrology_core.cpp

@@ -76 +76 @@
                                 substep=0;
                                 femmodel->parameters->SetParam(subdt,TimesteppingTimeStepEnum);
+                                femmodel->parameters->SetParam(subtime,TimeEnum);
 
                                 /*intermiedaries to deal with averaging*/
@@ -119 +120 @@
                                 /*averaging the stack*/
                                 femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,hydro_averaging);
-
                                 /*And reseting to global time*/
                                 femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);

issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp

@@ -9 +9 @@
 void InputUpdateFromSolutionx(FemModel* femmodel,Vector<IssmDouble>* solution){
 
-        /*Display message*/
-        if(VerboseModule()) _printf0_("   Updating inputs from solution\n");
-
         /*GetAnalysis*/
         int analysisenum;
         femmodel->parameters->FindParam(&analysisenum,AnalysisTypeEnum);
         Analysis* analysis = EnumToAnalysis(analysisenum);
+
+        /*Display message*/
+        if(VerboseModule()) _printf0_("   Updating inputs from solution for " << EnumToStringx(analysisenum) << "\n");
 
         /*Get local vector with both masters and slaves:*/

issm/trunk-jpl/src/m/classes/hydrologydc.m

@@ +1 @@
+
 %Hydrologydc class definition
 %

issm/trunk-jpl/test/NightlyRun/test905.m

@@ -33 +33 @@
 
 md.hydrology.steps_per_step=5;
-md.smb.steps_per_step=7.;
-md.timestepping.time_step=2.;
+md.smb.steps_per_step=10;
+md.timestepping.time_step=1.;
 md.timestepping.final_time=20.0;
 
@@ -59 +59 @@
              'SedimentWaterHead5','EplWaterHead5','SedimentHeadResidual5',...
              'SedimentWaterHead9','EplWaterHead9','SedimentHeadResidual9',...
-             'EplWaterHead10','EplWaterHeadSubstep10','SedimentWaterHead10',...
-             'SedimentWaterHeadSubstep10'}
+             'EplWaterHead20','EplWaterHeadSubstep20','SedimentWaterHead20',...
+             'SedimentWaterHeadSubstep20'}
 field_tolerances={1e-13,1e-13,1e-13,...
                   1e-13,1e-13,1e-13,...

issm/trunk-jpl/test/NightlyRun/test905.py

@@ -42 +42 @@
 md.initialization.epl_thickness = np.ones((md.mesh.numberofvertices))
 
+#try to keep the different steps as common multipliers of eachother
+#for the sake of precision the values of the model used as input should be on a shorter time-step-
+#you can plot the results of this test and check how the time stepping is doing (see bellow)
 md.hydrology.steps_per_step = 5
-md.smb.steps_per_step = 7  #md.hydrology.steps_per_step
-md.timestepping.time_step = 2.
+md.smb.steps_per_step = 10  #md.hydrology.steps_per_step
+md.timestepping.time_step = 1.
 md.timestepping.final_time = 20.0
 
 smb_step = md.timestepping.time_step / md.smb.steps_per_step
+hydro_step = md.timestepping.time_step / md.hydrology.steps_per_step
 duration = np.arange(md.timestepping.start_time, md.timestepping.final_time + smb_step, smb_step)
 
@@ -55 +59 @@
 md.smb.accugrad = 0.0
 
+#md.smb.runoffref = 9. * ddf * np.ones(np.shape(duration))  #constant input for testing
 md.smb.runoffref = 0.9 * duration * ddf
 md.smb.runoffref = np.vstack((md.smb.runoffref, duration))
@@ -64 +69 @@
 
 md = solve(md, 'Transient')
-
 field_names = ['SedimentWaterHead1', 'EplWaterHead1', 'SedimentHeadResidual1',
                'SedimentWaterHead4', 'EplWaterHead4', 'SedimentHeadResidual4',
                'SedimentWaterHead5', 'EplWaterHead5', 'SedimentHeadResidual5',
                'SedimentWaterHead9', 'EplWaterHead9', 'SedimentHeadResidual9',
-               'EplWaterHead10', 'EplWaterHeadSubstep10', 'SedimentWaterHead10',
-               'SedimentWaterHeadSubstep10']
+               'EplWaterHead20', 'EplWaterHeadSubstep20', 'SedimentWaterHead20',
+               'SedimentWaterHeadSubstep20']
 field_tolerances = [1e-13, 1e-13, 1e-13,
                     1e-13, 1e-13, 1e-13,
@@ -93 +97 @@
                 md.results.TransientSolution[-1].SedimentHead,
                 md.results.TransientSolution[-1].SedimentHeadSubstep]
+
+
+#===Stepping and plotting reference
+# import matplotlib as mpl
+# import matplotlib.pyplot as plt
+
+# agregatedInput = np.zeros(len(md.results.TransientSolution))
+# agregatedInputSub = np.zeros(len(md.results.TransientSolution))
+# RefInput = np.zeros(len(md.results.TransientSolution))
+# RefInputSub = np.zeros(len(md.results.TransientSolution))
+# Input = np.zeros(len(md.results.TransientSolution))
+# InputSub = np.zeros(len(md.results.TransientSolution))
+# SedVol = np.zeros(len(md.results.TransientSolution))
+# SedVolSub = np.zeros(len(md.results.TransientSolution))
+# EplVol = np.zeros(len(md.results.TransientSolution))
+# EplVolSub = np.zeros(len(md.results.TransientSolution))
+# SmbTimer = np.zeros(len(md.results.TransientSolution))
+# HydroTimer = np.zeros(len(md.results.TransientSolution))
+# TimerSub = np.zeros(len(md.results.TransientSolution))
+# HeadEVol = np.zeros(len(md.results.TransientSolution))
+# sedstore = md.materials.rho_freshwater * md.constants.g * md.hydrology.sediment_porosity * md.hydrology.sediment_thickness * ((md.hydrology.sediment_compressibility / md.hydrology.sediment_porosity) + md.hydrology.water_compressibility)
+# eplstore = md.materials.rho_freshwater * md.constants.g * md.hydrology.sediment_porosity * ((md.hydrology.sediment_compressibility / md.hydrology.sediment_porosity) + md.hydrology.water_compressibility)
+
+# for i, res in enumerate(md.results.TransientSolution):
+
+#     TimerSub[i] = res.time
+#     RefInputSub[i] = max(0, res.time * 0.9 * ddf + (-6.5e-3 * 1000 * ddf))
+#     try:
+#         InputSub[i] = np.nanmean(res.SmbRunoffSubstep)
+#         SedVolSub[i] = np.nanmean(res.SedimentHeadSubstep) * sedstore
+#         EplVolSub[i] = np.nanmean(res.EplHeadSubstep) * eplstore * np.nanmean(res.HydrologydcEplThicknessSubstep)
+#         substep = True
+#         if i > 0:
+#             agregatedInputSub[i:] += 0.5 * ((1 + smb_step) * np.nanmean(res.SmbRunoffSubstep) + (1 - smb_step) * np.nanmean(md.results.TransientSolution[i - 1].SmbRunoffSubstep)) * (TimerSub[i] - TimerSub[i - 1])
+#         else:
+#             agregatedInputSub[i:] += np.nanmean(res.SmbRunoffSubstep) * TimerSub[i]
+
+#     except AttributeError:
+#         substep = False
+#     SmbTimer[i] = res.time - (0.5 * md.timestepping.time_step * (1 - smb_step))
+#     HydroTimer[i] = res.time - (0.5 * md.timestepping.time_step * (1 - hydro_step))
+
+#     RefInput[i] = max(0, (SmbTimer[i] * 0.9 * ddf + (-6.5e-3 * 1000 * ddf)))
+#     Input[i] = np.nanmean(res.SmbRunoff)
+#     SedVol[i] = np.nanmean(res.SedimentHead) * sedstore
+#     EplVol[i] = np.nanmean(res.EplHead) * eplstore * np.nanmean(res.HydrologydcEplThickness)
+#     if i > 0:
+#         agregatedInput[i:] += 0.5 * ((1 + hydro_step) * np.nanmean(res.SmbRunoff) + (1 - hydro_step) * np.nanmean(md.results.TransientSolution[i - 1].SmbRunoff)) * (SmbTimer[i] - SmbTimer[i - 1])
+#     else:
+#         agregatedInput[i:] += np.nanmean(res.SmbRunoff) * SmbTimer[i]
+
+
+# layout, ax = plt.subplots(2, 2, sharex=True, sharey=False, figsize=(10, 10))
+
+# DiffTimer = (TimerSub[:-1] + 0.5 * (smb_step + hydro_step))
+# Indiff = (DiffTimer) * 0.9 * ddf + (-6.5e-3 * 1000 * ddf)
+# Indiff[np.where(Indiff < 0)] = 0
+
+# ax[0, 0].plot(SmbTimer, RefInput - Input, 'g+')
+# ax[0, 0].plot(DiffTimer, (Indiff - np.diff((SedVol+EplVol)) / md.timestepping.time_step), 'kx')
+
+# ax[0, 1].plot(SmbTimer, RefInput, 'm+')
+# ax[0, 1].plot(SmbTimer, Input, 'g+')
+# ax[0, 1].plot(DiffTimer, np.diff(SedVol+EplVol) / md.timestepping.time_step, 'r')
+
+# ax[1, 0].plot(HydroTimer, (agregatedInput - SedVol - EplVol), 'b+')
+
+# ax[1, 1].plot(HydroTimer, SedVol + EplVol, 'r+')
+# ax[1, 1].plot(SmbTimer, agregatedInput, 'b+')
+
+# if substep:
+#     SubDiffTimer = (SmbTimer[1:])
+#     SubIndiff = (SubDiffTimer) * 0.9 * ddf + (-6.5e-3 * 1000 * ddf)
+#     SubIndiff[np.where(SubIndiff < 0)] = 0
+
+#     ax[0, 0].plot(TimerSub, RefInputSub - InputSub, 'm')
+#     ax[0, 0].plot(SubDiffTimer, (SubIndiff - np.diff(SedVolSub + EplVolSub) / md.timestepping.time_step), 'y')
+
+#     ax[0, 1].plot(TimerSub, RefInputSub, 'm')
+#     ax[0, 1].plot(TimerSub, InputSub, 'g')
+
+#     ax[1, 0].plot(TimerSub, (agregatedInputSub - SedVolSub - EplVolSub), 'r')
+
+#     ax[1, 1].plot(TimerSub, SedVolSub + EplVolSub, 'r')
+#     ax[1, 1].plot(TimerSub, SedVolSub, 'r', alpha=0.5)
+#     ax[1, 1].plot(TimerSub, EplVolSub, 'r', alpha=0.5)
+#     ax[1, 1].plot(TimerSub, agregatedInputSub, 'b')