Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 25023) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 25024) @@ -299,4 +299,7 @@ /*Intermediaries */ int smb_model; + int smb_averaging; + int smbsubstepping; + int hydrologysubstepping; IssmDouble dt,scalar,water_head; IssmDouble water_load,transfer,runoff_value; @@ -307,4 +310,5 @@ IssmDouble *xyz_list = NULL; Input2 *old_wh_input = NULL; + Input2 *dummy_input = NULL; Input2 *surface_runoff_input = NULL; @@ -321,4 +325,5 @@ 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,9 +334,27 @@ 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(dummy_input); _assert_(surface_runoff_input); } Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 25023) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 25024) @@ -375,5 +375,4 @@ basalelement->FindParam(&time,TimeEnum); - if(dt!= 0.){ old_wh_input = basalelement->GetInput2(SedimentHeadOldEnum); _assert_(old_wh_input); @@ -384,10 +383,13 @@ 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); } Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 25023) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 25024) @@ -1009,133 +1009,13 @@ /*}}}*/ 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(×teps,&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_timeGetPentaInput(0); - - int interpolation = input->GetInterpolation(); - _assert_(interpolation==P1Enum); - /*Intermediaries*/ - int numindices; - int indices[NUMVERTICES]; - numindices = NUMVERTICES; - for(int i=0;ilid; - input->Serve(numindices,&indices[0]); - - for(int iv=0;ivGaussVertex(iv); - input->GetInputValue(¤t_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;ilid; - input->Serve(numindices,&indices[0]); - - for(int iv=0;ivGaussVertex(iv); - input->GetInputValue(¤t_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;iv1.e-50); - averaged_values[iv] = dt*1./current_values[iv]; - } - } - else{ - for(int iv=0;iv1.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;ivAddInput2(averagedinput_enum,&averaged_values[0],P1Enum); - - /*Cleanup*/ - delete gauss; - xDelete(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); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25023) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25024) @@ -914,5 +914,5 @@ } /*}}}*/ -void Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/ +void Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/ IssmDouble xyz_list[NUMVERTICES][3]; @@ -1918,5 +1918,9 @@ /*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,113 +2118,11 @@ _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(×teps,&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_timeGetTriaInput(0); - _assert_(input->GetInterpolation()==P1Enum); - this->InputServe(input); - for(int iv=0;ivGaussVertex(iv); - input->GetInputValue(¤t_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;ivGaussVertex(iv); - input->GetInputValue(¤t_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;iv1.e-50); - averaged_values[iv] = dt*1./current_values[iv]; - } - } - else{ - for(int iv=0;iv1.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;ivAddInput2(averagedinput_enum,&averaged_values[0],P1Enum); - - /*Cleanup*/ - delete gauss; - xDelete(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,5 +5231,5 @@ rho_ice=FindParam(MaterialsRhoIceEnum); rho_earth=FindParam(MaterialsEarthDensityEnum); - + /*recover earth area: */ this->parameters->FindParam(&planetarea,SealevelPlanetAreaEnum); @@ -5356,5 +5258,5 @@ } - // 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,5 +5264,5 @@ 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,5 +5471,5 @@ 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,5 +5497,5 @@ 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,5 +5518,5 @@ IssmDouble* indices=NULL; #endif - + /*elastic green function:*/ int index; @@ -5656,5 +5558,5 @@ GE=xNewZeroInit(gsize); } - + /*compute area:*/ area=GetAreaSpherical(); @@ -5676,5 +5578,5 @@ constant=3/rho_earth*area/planetarea; - + for(int i=0;i(M-1); index=reCast(indices[i]); - + /*Rigid earth gravitational perturbation: */ if(computerigid){ @@ -5695,5 +5597,5 @@ G[i]+=G_elastic_precomputed[index]; } - G[i]=G[i]*constant; + G[i]=G[i]*constant; /*Elastic components:*/ @@ -5728,10 +5630,10 @@ /*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,5 +5654,5 @@ int bp_compute_fingerprints= 0; - + /*Compute bottom pressure contribution from ocean if requested:*/ this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum); @@ -5795,8 +5697,8 @@ /*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,7 +5707,7 @@ /*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,5 +5722,5 @@ /*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,10 +5728,10 @@ /*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,5 +5794,5 @@ /*elastic green function:*/ IssmDouble* G=NULL; - + /*water properties: */ IssmDouble rho_water; @@ -5908,5 +5810,5 @@ /*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,12 +5845,12 @@ /*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,5 +5858,5 @@ /*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,8 +5906,8 @@ /*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); } Index: /issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp (revision 25023) +++ /issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp (revision 25024) @@ -388,4 +388,20 @@ } }/*}}}*/ +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(input)->GetPentaInput(start_time,end_time,averaging_method); + } + else{ + _error_("Input "<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(this->current_input); + +} +/*}}}*/ + void TransientInput2::SetCurrentTimeInput(IssmDouble time){/*{{{*/ @@ -404,5 +419,6 @@ /*If already processed return*/ - if(this->current_step>this_step-1.e-5 && this->current_stepcurrent_step-this_step)<1.e-5) return; + // if(this->current_step>this_step-1.e-5 && this->current_steptimesteps,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(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_timetimesteps[0]); - } - else if(offset==this->numtimesteps-1)dt=end_time-this->timesteps[offset]; - else if(offset==start_offset && this->timesteps[offset]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,5 +487,5 @@ case 0: /*Arithmetic mean*/ if(offset==start_offset){ - this->current_input = stepinput->copy(); + this->current_input=stepinput->copy(); this->current_input->Scale(dt); } Index: /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h (revision 25023) +++ /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h (revision 25024) @@ -55,4 +55,5 @@ 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); Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 25023) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 25024) @@ -76,4 +76,5 @@ substep=0; femmodel->parameters->SetParam(subdt,TimesteppingTimeStepEnum); + femmodel->parameters->SetParam(subtime,TimeEnum); /*intermiedaries to deal with averaging*/ @@ -119,5 +120,4 @@ /*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); Index: /issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp (revision 25023) +++ /issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp (revision 25024) @@ -9,11 +9,11 @@ void InputUpdateFromSolutionx(FemModel* femmodel,Vector* 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:*/ Index: /issm/trunk-jpl/src/m/classes/hydrologydc.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologydc.m (revision 25023) +++ /issm/trunk-jpl/src/m/classes/hydrologydc.m (revision 25024) @@ -1,2 +1,3 @@ + %Hydrologydc class definition % Index: /issm/trunk-jpl/test/NightlyRun/test905.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test905.m (revision 25023) +++ /issm/trunk-jpl/test/NightlyRun/test905.m (revision 25024) @@ -33,6 +33,6 @@ 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,6 +59,6 @@ '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,... Index: /issm/trunk-jpl/test/NightlyRun/test905.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test905.py (revision 25023) +++ /issm/trunk-jpl/test/NightlyRun/test905.py (revision 25024) @@ -42,10 +42,14 @@ 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,4 +59,5 @@ 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,11 +69,10 @@ 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,2 +97,89 @@ 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')