Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 27176) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 27177) @@ -11,8 +11,13 @@ void HydrologyDCEfficientAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + bool isefficientlayer; int hydrology_model; int eplflip_lock; int eplthickcomp; - bool isefficientlayer; + IssmDouble eplinitthick; + IssmDouble eplcolapsethick; + IssmDouble eplmaxthick; + IssmDouble eplcond; + /*retrieve some parameters: */ iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); @@ -27,10 +32,20 @@ if(!isefficientlayer) return; - /*If yes, initialize a flip flop counter*/ + /*If yes, get the parameters*/ iomodel->FetchData(&eplflip_lock,"md.hydrology.eplflip_lock"); + iomodel->FetchData(&eplthickcomp,"md.hydrology.epl_thick_comp"); + parameters->AddObject(new IntParam(HydrologydcEplflipLockEnum,eplflip_lock)); - - iomodel->FetchData(&eplthickcomp,"md.hydrology.epl_thick_comp"); parameters->AddObject(new IntParam(HydrologydcEplThickCompEnum,eplthickcomp)); + + iomodel->FetchData(&eplinitthick,"md.hydrology.epl_initial_thickness"); + iomodel->FetchData(&eplcolapsethick,"md.hydrology.epl_colapse_thickness"); + iomodel->FetchData(&eplmaxthick,"md.hydrology.epl_max_thickness"); + iomodel->FetchData(&eplcond,"md.hydrology.epl_conductivity"); + parameters->AddObject(new DoubleParam(HydrologydcEplInitialThicknessEnum,eplinitthick)); + parameters->AddObject(new DoubleParam(HydrologydcEplColapseThicknessEnum,eplcolapsethick)); + parameters->AddObject(new DoubleParam(HydrologydcEplMaxThicknessEnum,eplmaxthick)); + parameters->AddObject(new DoubleParam(HydrologydcEplConductivityEnum,eplcond)); + }/*}}}*/ void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 27176) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 27177) @@ -23,8 +23,11 @@ int numoutputs; bool isefficientlayer; + bool sliceadapt; IssmDouble penalty_factor; IssmDouble rel_tol; IssmDouble leakagefactor; IssmDouble sedimentlimit; + IssmDouble sed_poro; + IssmDouble sed_thick; char** requestedoutputs = NULL; @@ -43,4 +46,5 @@ iomodel->FetchData(&hydro_maxiter, "md.hydrology.max_iter" ); iomodel->FetchData(&hydroslices, "md.hydrology.steps_per_step"); + iomodel->FetchData(&sliceadapt, "md.hydrology.step_adapt"); iomodel->FetchData(&averaging_method, "md.hydrology.averaging"); iomodel->FetchData(&isefficientlayer, "md.hydrology.isefficientlayer"); @@ -55,9 +59,16 @@ parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter)); parameters->AddObject(new IntParam(HydrologyStepsPerStepEnum,hydroslices)); + parameters->AddObject(new BoolParam(HydrologyStepAdaptEnum,sliceadapt)); parameters->AddObject(new IntParam(HydrologyAveragingEnum,averaging_method)); - parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer)); parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor)); parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol)); + + iomodel->FetchData(&sed_poro, "md.hydrology.sediment_porosity" ); + iomodel->FetchData(&sed_thick, "md.hydrology.sediment_thickness" ); + + parameters->AddObject(new DoubleParam(HydrologydcSedimentPorosityEnum,sed_poro)); + parameters->AddObject(new DoubleParam(HydrologydcSedimentThicknessEnum,sed_thick)); + if(transfer_flag==1){ iomodel->FetchData(&leakagefactor,"md.hydrology.leakage_factor"); @@ -110,4 +121,6 @@ iomodel->FetchDataToInput(inputs,elements,"md.hydrology.sediment_transmitivity",HydrologydcSedimentTransmitivityEnum); iomodel->FetchDataToInput(inputs,elements,"md.hydrology.mask_thawed_node",HydrologydcMaskThawedNodeEnum); + + if(iomodel->domaintype!=Domain2DhorizontalEnum){ iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 27176) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 27177) @@ -57,9 +57,9 @@ InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum); - /*grab tws from the hydrology.spcwatercolumn field input and update + /*grab tws from the hydrology.spcwatercolumn field input and update * the solution with it:*/ Vector* ug = NULL; GetVectorFromInputsx(&ug,femmodel,HydrologyTwsSpcEnum,VertexPIdEnum); - InputUpdateFromSolutionx(femmodel,ug); + InputUpdateFromSolutionx(femmodel,ug); /*solid earth considerations:*/ @@ -83,13 +83,17 @@ if(ThawedNodes>0){ /*check if we need sub steps*/ - int dtslices; - femmodel->parameters->FindParam(&dtslices,HydrologyStepsPerStepEnum); - - if(dtslices>1){ - int substep, numaveragedinput, hydro_averaging; - IssmDouble global_time, subtime, yts; + int dtslices; + bool sliceadapt; + bool conv_fail=false; + femmodel->parameters->FindParam(&dtslices,HydrologyStepsPerStepEnum); + femmodel->parameters->FindParam(&sliceadapt,HydrologyStepAdaptEnum); + + if(dtslices>1 || sliceadapt){ + int step, substep, numaveragedinput, hydro_averaging, remainingslices; + IssmDouble global_time, subtime, yts, remainingtime, averagetime; IssmDouble dt, subdt; femmodel->parameters->FindParam(&global_time,TimeEnum); + femmodel->parameters->FindParam(&step,StepEnum); femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); @@ -124,7 +128,9 @@ } femmodel->InitTransientInputx(&transientinput[0],numaveragedinput); + averagetime=0; while(substepparameters->SetParam(subtime,TimeEnum); @@ -138,7 +144,32 @@ } /*Proceed now to heads computations*/ - solutionsequence_hydro_nonlinear(femmodel); - /*If we have a sub-timestep we store the substep inputs in a transient input here*/ - femmodel->StackTransientInputonBasex(&substepinput[0],&transientinput[0],subtime,numaveragedinput); + solutionsequence_hydro_nonlinear(femmodel, &conv_fail); + if(conv_fail){ + /*convergence failed, we want to go back to the begining of the main step and increase the number of subslices*/ + /*First we get teh time and step counter back to the begining of the step that did not converge*/ + averagetime-=subtime*subdt; + subtime-=subdt; + substep-=1; + /*compute the number of slice that are remaining and the time left in the timestep*/ + remainingslices=dtslices-substep; + remainingtime=global_time-subtime; + /*We double the number of remaining slices and compute their duration*/ + dtslices=dtslices-remainingslices+(2*remainingslices); + subdt=remainingtime/(2*remainingslices); + _printf0_("convergence failed for sub-step "<< substep <<" total number of slice is now "<< dtslices <<" for step "<500){ + _error_(" We reached (" << dtslices << ") which exceeds the hard limit of 500"); + } + + femmodel->parameters->SetParam(subdt,TimesteppingTimeStepEnum); + femmodel->parameters->SetParam(subtime,TimeEnum); + + } + else{ + /*If we have a sub-timestep we store the substep inputs in a transient input here*/ + femmodel->StackTransientInputonBasex(&substepinput[0],&transientinput[0],subtime,numaveragedinput); + } } /*averaging the stack*/ @@ -147,4 +178,8 @@ femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum); femmodel->parameters->SetParam(global_time,TimeEnum); + if(save_results){ + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,HydrologySubstepsEnum,dtslices,step,global_time)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,HydrologySubTimeEnum,(averagetime/dt)/yts,step,global_time)); + } } else{ @@ -156,5 +191,5 @@ /*Proceed now to heads computations*/ if(VerboseSolution()) _printf0_(" computing water heads\n"); - solutionsequence_hydro_nonlinear(femmodel); + solutionsequence_hydro_nonlinear(femmodel, &conv_fail); /*If no substeps are present we want to duplicate the results for coupling purposes*/ InputDuplicatex(femmodel,SedimentHeadSubstepEnum,SedimentHeadEnum); Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 27176) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 27177) @@ -10,5 +10,5 @@ #include "../solutionsequences/solutionsequences.h" -void solutionsequence_hydro_nonlinear(FemModel* femmodel){ +void solutionsequence_hydro_nonlinear(FemModel* femmodel, bool* pconv_fail){ /*solution : */ Vector* ug_sed=NULL; @@ -16,14 +16,14 @@ Vector* uf_sed_sub_iter=NULL; Vector* ug_sed_main_iter=NULL; + Vector* ug_sed_init=NULL; Vector* ug_epl=NULL; Vector* uf_epl=NULL; Vector* uf_epl_sub_iter=NULL; - Vector* ug_epl_sub_iter=NULL; Vector* ug_epl_main_iter=NULL; + Vector* ug_epl_init=NULL; Vector* ys=NULL; Vector* dug=NULL; - Vector* duf=NULL; Matrix* Kff=NULL; @@ -38,10 +38,12 @@ bool sedconverged,eplconverged,hydroconverged; bool isefficientlayer; + bool sliceadapt; int constraints_converged; int num_unstable_constraints; int sedcount,eplcount,hydrocount; int hydro_maxiter; + int epl_fsize,epl_sub_fsize,epl_main_fsize; IssmDouble sediment_kmax; - IssmDouble eps_hyd; + IssmDouble eps_res,eps_rel,eps_abs; IssmDouble ndu_sed,nu_sed; IssmDouble ndu_epl,nu_epl; @@ -50,6 +52,9 @@ femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + femmodel->parameters->FindParam(&sliceadapt,HydrologyStepAdaptEnum); femmodel->parameters->FindParam(&hydro_maxiter,HydrologydcMaxIterEnum); - femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum); + femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum); + femmodel->parameters->FindParam(&eps_rel,HydrologydcRelTolEnum); + femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum); hydrocount=1; hydroconverged=false; @@ -59,5 +64,11 @@ /*{{{*//*Retrieve inputs as the initial state for the non linear iteration*/ GetBasalSolutionFromInputsx(&ug_sed,femmodel); + /*Initialize the IDS element mask to exclude frozen nodes*/ + inefanalysis->ElementizeIdsMask(femmodel); + Reducevectorgtofx(&uf_sed, ug_sed, femmodel->nodes,femmodel->parameters); + ug_sed_init=ug_sed->Duplicate(); + ug_sed->Copy(ug_sed_init); + if(isefficientlayer) { inefanalysis = new HydrologyDCInefficientAnalysis(); @@ -65,15 +76,13 @@ femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); GetBasalSolutionFromInputsx(&ug_epl,femmodel); - /*Initialize the EPL element mask*/ inefanalysis->ElementizeEplMask(femmodel); effanalysis->InitZigZagCounter(femmodel); - /*Initialize the IDS element mask*/ - femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); - inefanalysis->ElementizeIdsMask(femmodel); + Reducevectorgtofx(&uf_epl, ug_epl, femmodel->nodes,femmodel->parameters); + ug_epl_init=ug_epl->Duplicate(); + ug_epl->Copy(ug_epl_init); } /*}}}*/ /*The real computation starts here, outermost loop is on the two layer system*/ for(;;){ - sedcount=1; eplcount=1; @@ -105,4 +114,5 @@ /*{{{*//*Core of the computation*/ if(VerboseSolution()) _printf0_("Building Sediment Matrix...\n"); + femmodel->profiler->Start(SEDMatrix); SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax,femmodel); @@ -114,7 +124,6 @@ femmodel->profiler->Start(SOLVER); Solverx(&uf_sed,Kff,pf,uf_sed_sub_iter,df,femmodel->parameters); + delete df; femmodel->profiler->Stop(SOLVER); - - delete Kff; delete pf; delete df; delete ug_sed; femmodel->profiler->Start(SEDUpdate); @@ -125,8 +134,18 @@ /*}}}*/ if (!sedconverged){ + /*First check that all the penalizations are applied*/ if(VerboseConvergence()) _printf0_(" # Sediment unstable constraints = " << num_unstable_constraints << "\n"); - if(num_unstable_constraints==0) sedconverged = true; + if(num_unstable_constraints==0) { + sedconverged = true; + } + else{//clean up + delete Kff; + delete pf; + } if (sedcount>=hydro_maxiter){ + delete ug_sed;delete uf_sed;delete inefanalysis; delete ug_sed_main_iter; + if(isefficientlayer)delete ug_epl;delete uf_epl;delete effanalysis; delete ug_epl_main_iter; _error_(" maximum number of Sediment iterations (" << hydro_maxiter << ") exceeded"); + } } @@ -135,23 +154,9 @@ if(sedconverged)break; } - /*}}}*//*End of the sediment penalization loop*/ sedconverged=false; - /*Checking convergence on the value of the sediment head*/ - duf=uf_sed_sub_iter->Duplicate();_assert_(duf); - uf_sed_sub_iter->Copy(duf); - duf->AYPX(uf_sed,-1.0); - ndu_sed=duf->Norm(NORM_TWO); - delete duf; - nu_sed=uf_sed_sub_iter->Norm(NORM_TWO); - if (xIsNan(ndu_sed) || xIsNan(nu_sed)) _error_("convergence criterion is NaN!"); - if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the layer is empty*/ - if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << "%, aiming lower than " << eps_hyd*10*100 << " %\n"); - if((ndu_sed/nu_sed)parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum); @@ -176,6 +181,4 @@ for(;;){ eplconverged=false; - ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter); - ug_epl->Copy(ug_epl_sub_iter); /*{{{*//*Retrieve the EPL head slopes and compute EPL Thickness*/ if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); @@ -195,19 +198,17 @@ femmodel->profiler->Stop(EPLMasking); if(VerboseSolution()) _printf0_("Building EPL Matrix...\n"); + uf_epl_sub_iter=uf_epl->Duplicate();_assert_(uf_epl_sub_iter); + uf_epl->Copy(uf_epl_sub_iter); + uf_epl->GetSize(&epl_sub_fsize); femmodel->profiler->Start(EPLMatrices); SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); CreateNodalConstraintsx(&ys,femmodel->nodes); - Reduceloadx(pf,Kfs,ys); delete Kfs; - delete uf_epl; + Reduceloadx(pf,Kfs,ys); + delete Kfs;delete uf_epl; femmodel->profiler->Stop(EPLMatrices); femmodel->profiler->Start(SOLVER); Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters); femmodel->profiler->Stop(SOLVER); - - delete Kff; delete pf; delete df; - delete uf_epl_sub_iter; - uf_epl_sub_iter=uf_epl->Duplicate();_assert_(uf_epl_sub_iter); - uf_epl->Copy(uf_epl_sub_iter); - delete ug_epl; + delete df;delete ug_epl; femmodel->profiler->Start(EPLUpdate); Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys; @@ -215,25 +216,34 @@ ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); femmodel->profiler->Stop(EPLUpdate); - - dug=ug_epl_sub_iter->Duplicate();_assert_(dug); - ug_epl_sub_iter->Copy(dug); - dug->AYPX(ug_epl,-1.0); - ndu_epl=dug->Norm(NORM_TWO); - delete dug; - nu_epl=ug_epl_sub_iter->Norm(NORM_TWO); - if (xIsNan(ndu_epl) || xIsNan(nu_epl)) _error_("convergence criterion is NaN!"); - if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/ - if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner EPL Convergence criterion:" << ndu_epl/nu_epl*100 << "%, aiming lower than " << eps_hyd*10*100 << " %\n"); - if((ndu_epl/nu_epl)=hydro_maxiter){ + uf_epl->GetSize(&epl_fsize); + if(epl_fsize-epl_sub_fsize==0){ + convergence(&eplconverged,Kff,pf,uf_epl,uf_epl_sub_iter,eps_res,eps_rel,eps_abs); + delete Kff; delete pf; + /* if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /\*Hacking the case where the EPL is used but empty*\/ */ + } + else{ + delete Kff; delete pf; + } + + if (eplcount>=hydro_maxiter*9/10 && sliceadapt && !eplconverged) { + if(VerboseSolution()) _printf0_("epl did not converged after "<SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); + InputUpdateFromSolutionx(femmodel,ug_sed_init); + delete ug_sed_init; + break; + } + else if (eplcount>=hydro_maxiter){ + delete ug_sed;delete uf_sed;delete inefanalysis;delete ug_sed_main_iter; + delete ug_epl;delete uf_epl;delete effanalysis;delete ug_epl_main_iter; _error_(" maximum number of EPL iterations (" << hydro_maxiter << ") exceeded"); + } eplcount++; - - delete ug_epl_sub_iter; + delete uf_epl_sub_iter; if(eplconverged){ if(VerboseSolution()) _printf0_("eplconverged...\n"); - InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum); - InputUpdateFromSolutionx(femmodel,ug_epl); effanalysis->ResetCounter(femmodel); break; @@ -242,4 +252,7 @@ } femmodel->profiler->Stop(EPLLOOP); + if(*pconv_fail){ + break; + } /*}}}*//*End of the global EPL loop*/ /*{{{*//*Now dealing with the convergence of the whole system*/ @@ -264,18 +277,32 @@ if (xIsNan(ndu_epl) || xIsNan(nu_epl)) _error_("EPL convergence criterion is NaN!"); if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/ - if (!xIsNan(eps_hyd)){ - if ((ndu_epl/nu_epl)(eps_rel)){ + if ((ndu_epl/nu_epl)=hydro_maxiter*9/10 && sliceadapt && !hydroconverged) { + if(VerboseSolution()) _printf0_("hydrology main loop did not converged after "<SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); + InputUpdateFromSolutionx(femmodel,ug_sed_init); + delete ug_sed_init; + break; + } if (hydrocount>=hydro_maxiter){ _error_(" maximum number for hydrological global iterations (" << hydro_maxiter << ") exceeded"); + delete ug_sed;delete uf_sed;delete effanalysis; + delete ug_epl; delete uf_epl; delete inefanalysis; } } @@ -284,14 +311,14 @@ } /*}}}*/ - if(isefficientlayer)InputUpdateFromSolutionx(femmodel,ug_epl); - femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); - InputUpdateFromSolutionx(femmodel,ug_sed); + /*To deal with adaptative stepping we only save results if we are actually converged*/ + if(hydroconverged){ + if(isefficientlayer){ + delete ug_epl_init; + } + delete ug_sed_init; + } /*Free resources: */ - delete ug_epl; - delete ug_sed; - delete uf_sed; - delete uf_epl; - delete uf_epl_sub_iter; - delete inefanalysis; - delete effanalysis; + delete ug_epl;delete ug_sed; + delete uf_sed; delete uf_epl; + delete inefanalysis; delete effanalysis; } Index: /issm/trunk-jpl/src/m/classes/hydrologydc.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologydc.m (revision 27176) +++ /issm/trunk-jpl/src/m/classes/hydrologydc.m (revision 27177) @@ -14,4 +14,5 @@ max_iter = 0; steps_per_step = 0; + step_adapt = 0; averaging = 0; sedimentlimit_flag = 0; @@ -66,5 +67,5 @@ list=[list,{'EplHead','HydrologydcMaskEplactiveNode','HydrologydcMaskEplactiveElt','EplHeadSlopeX','EplHeadSlopeY','HydrologydcEplThickness'}]; end - if self.steps_per_step>1, + if self.steps_per_step>1 | self.step_adapt, list = [list,'EffectivePressureSubstep','SedimentHeadSubstep']; if self.isefficientlayer, @@ -90,4 +91,5 @@ self.max_iter = 100; self.steps_per_step = 1; + self.step_adapt = 0; self.averaging = 0; self.sedimentlimit_flag = 0; @@ -126,4 +128,5 @@ md = checkfield(md,'fieldname','hydrology.max_iter','>',0,'numel',1); md = checkfield(md,'fieldname','hydrology.steps_per_step','>',0,'numel',1); + md = checkfield(md,'fieldname','hydrology.step_adapt','numel',1,'values',[0 1]); md = checkfield(md,'fieldname','hydrology.averaging','numel',[1],'values',[0 1 2]); md = checkfield(md,'fieldname','hydrology.sedimentlimit_flag','numel',[1],'values',[0 1 2 3]); @@ -172,4 +175,5 @@ fielddisplay(self,'max_iter','maximum number of nonlinear iteration'); fielddisplay(self,'steps_per_step','number of hydrology steps per timestep'); + fielddisplay(self,'step_adapt', 'adaptative sub stepping [1: true 0: false] default is 0'); fielddisplay(self, 'averaging', 'averaging methods from short to long steps'); disp(sprintf('%55s 0: Arithmetic (default)')); @@ -226,4 +230,5 @@ WriteData(fid,prefix,'object',self,'fieldname','max_iter','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','step_adapt','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','sedimentlimit_flag','format','Integer'); Index: /issm/trunk-jpl/src/m/classes/hydrologydc.py =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologydc.py (revision 27176) +++ /issm/trunk-jpl/src/m/classes/hydrologydc.py (revision 27177) @@ -21,4 +21,5 @@ self.max_iter = 0 self.steps_per_step = 0 + self.step_adapt = 0 self.averaging = 0 self.sedimentlimit_flag = 0 @@ -53,5 +54,5 @@ def __repr__(self): # {{{ # TODO: - # - Convert all formatting to calls to .format (see any + # - Convert all formatting to calls to .format (see any # already converted .__repr__ method for examples) # @@ -65,4 +66,5 @@ string = "%s\n%s" % (string, fielddisplay(self, 'max_iter', 'maximum number of nonlinear iteration')) string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of hydrology steps per time step')) + string = "%s\n%s" % (string, fielddisplay(self, 'step_adapt', 'adaptative sub stepping [1: true 0: false] default is 0')) string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps')) string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)') @@ -134,4 +136,5 @@ self.max_iter = 100 self.steps_per_step = 1 + self.step_adapt = 0 self.averaging = 0 self.sedimentlimit_flag = 0 @@ -162,5 +165,5 @@ if self.isefficientlayer == 1: list.extend(['EplHead', 'HydrologydcMaskEplactiveNode', 'HydrologydcMaskEplactiveElt', 'EplHeadSlopeX', 'EplHeadSlopeY', 'HydrologydcEplThickness']) - if self.steps_per_step > 1: + if self.steps_per_step > 1 or self.step_adapt: list.extend(['EffectivePressureSubstep', 'SedimentHeadSubstep']) if self.isefficientlayer == 1: @@ -190,4 +193,5 @@ md = checkfield(md, 'fieldname', 'hydrology.max_iter', '>', 0., 'numel', [1]) md = checkfield(md, 'fieldname', 'hydrology.steps_per_step', '>=', 1, 'numel', [1]) + md = checkfield(md, 'fieldname', 'hydrology.step_adapt', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'hydrology.averaging', 'numel', [1], 'values', [0, 1, 2]) md = checkfield(md, 'fieldname', 'hydrology.sedimentlimit_flag', 'numel', [1], 'values', [0, 1, 2, 3]) @@ -233,4 +237,5 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'max_iter', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'step_adapt', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'sedimentlimit_flag', 'format', 'Integer')