![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 27177
- Timestamp:
- 08/02/22 02:54:18 (3 years ago)
- Location:
- issm/trunk-jpl/src
- Files:
-
- 6 edited
-
c/analyses/HydrologyDCEfficientAnalysis.cpp (modified) (2 diffs)
-
c/analyses/HydrologyDCInefficientAnalysis.cpp (modified) (4 diffs)
-
c/cores/hydrology_core.cpp (modified) (6 diffs)
-
c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (modified) (16 diffs)
-
m/classes/hydrologydc.m (modified) (6 diffs)
-
m/classes/hydrologydc.py (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
r27102 r27177 11 11 void HydrologyDCEfficientAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ 12 12 13 bool isefficientlayer; 13 14 int hydrology_model; 14 15 int eplflip_lock; 15 16 int eplthickcomp; 16 bool isefficientlayer; 17 IssmDouble eplinitthick; 18 IssmDouble eplcolapsethick; 19 IssmDouble eplmaxthick; 20 IssmDouble eplcond; 21 17 22 /*retrieve some parameters: */ 18 23 iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); … … 27 32 if(!isefficientlayer) return; 28 33 29 /*If yes, initialize a flip flop counter*/34 /*If yes, get the parameters*/ 30 35 iomodel->FetchData(&eplflip_lock,"md.hydrology.eplflip_lock"); 36 iomodel->FetchData(&eplthickcomp,"md.hydrology.epl_thick_comp"); 37 31 38 parameters->AddObject(new IntParam(HydrologydcEplflipLockEnum,eplflip_lock)); 32 33 iomodel->FetchData(&eplthickcomp,"md.hydrology.epl_thick_comp");34 39 parameters->AddObject(new IntParam(HydrologydcEplThickCompEnum,eplthickcomp)); 40 41 iomodel->FetchData(&eplinitthick,"md.hydrology.epl_initial_thickness"); 42 iomodel->FetchData(&eplcolapsethick,"md.hydrology.epl_colapse_thickness"); 43 iomodel->FetchData(&eplmaxthick,"md.hydrology.epl_max_thickness"); 44 iomodel->FetchData(&eplcond,"md.hydrology.epl_conductivity"); 45 parameters->AddObject(new DoubleParam(HydrologydcEplInitialThicknessEnum,eplinitthick)); 46 parameters->AddObject(new DoubleParam(HydrologydcEplColapseThicknessEnum,eplcolapsethick)); 47 parameters->AddObject(new DoubleParam(HydrologydcEplMaxThicknessEnum,eplmaxthick)); 48 parameters->AddObject(new DoubleParam(HydrologydcEplConductivityEnum,eplcond)); 49 35 50 }/*}}}*/ 36 51 void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ -
issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
r27102 r27177 23 23 int numoutputs; 24 24 bool isefficientlayer; 25 bool sliceadapt; 25 26 IssmDouble penalty_factor; 26 27 IssmDouble rel_tol; 27 28 IssmDouble leakagefactor; 28 29 IssmDouble sedimentlimit; 30 IssmDouble sed_poro; 31 IssmDouble sed_thick; 29 32 char** requestedoutputs = NULL; 30 33 … … 43 46 iomodel->FetchData(&hydro_maxiter, "md.hydrology.max_iter" ); 44 47 iomodel->FetchData(&hydroslices, "md.hydrology.steps_per_step"); 48 iomodel->FetchData(&sliceadapt, "md.hydrology.step_adapt"); 45 49 iomodel->FetchData(&averaging_method, "md.hydrology.averaging"); 46 50 iomodel->FetchData(&isefficientlayer, "md.hydrology.isefficientlayer"); … … 55 59 parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter)); 56 60 parameters->AddObject(new IntParam(HydrologyStepsPerStepEnum,hydroslices)); 61 parameters->AddObject(new BoolParam(HydrologyStepAdaptEnum,sliceadapt)); 57 62 parameters->AddObject(new IntParam(HydrologyAveragingEnum,averaging_method)); 58 59 63 parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer)); 60 64 parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor)); 61 65 parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol)); 66 67 iomodel->FetchData(&sed_poro, "md.hydrology.sediment_porosity" ); 68 iomodel->FetchData(&sed_thick, "md.hydrology.sediment_thickness" ); 69 70 parameters->AddObject(new DoubleParam(HydrologydcSedimentPorosityEnum,sed_poro)); 71 parameters->AddObject(new DoubleParam(HydrologydcSedimentThicknessEnum,sed_thick)); 72 62 73 if(transfer_flag==1){ 63 74 iomodel->FetchData(&leakagefactor,"md.hydrology.leakage_factor"); … … 110 121 iomodel->FetchDataToInput(inputs,elements,"md.hydrology.sediment_transmitivity",HydrologydcSedimentTransmitivityEnum); 111 122 iomodel->FetchDataToInput(inputs,elements,"md.hydrology.mask_thawed_node",HydrologydcMaskThawedNodeEnum); 123 124 112 125 if(iomodel->domaintype!=Domain2DhorizontalEnum){ 113 126 iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); -
issm/trunk-jpl/src/c/cores/hydrology_core.cpp
r27102 r27177 57 57 InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum); 58 58 59 /*grab tws from the hydrology.spcwatercolumn field input and update 59 /*grab tws from the hydrology.spcwatercolumn field input and update 60 60 * the solution with it:*/ 61 61 Vector<IssmDouble>* ug = NULL; 62 62 GetVectorFromInputsx(&ug,femmodel,HydrologyTwsSpcEnum,VertexPIdEnum); 63 InputUpdateFromSolutionx(femmodel,ug); 63 InputUpdateFromSolutionx(femmodel,ug); 64 64 65 65 /*solid earth considerations:*/ … … 83 83 if(ThawedNodes>0){ 84 84 /*check if we need sub steps*/ 85 int dtslices; 86 femmodel->parameters->FindParam(&dtslices,HydrologyStepsPerStepEnum); 87 88 if(dtslices>1){ 89 int substep, numaveragedinput, hydro_averaging; 90 IssmDouble global_time, subtime, yts; 85 int dtslices; 86 bool sliceadapt; 87 bool conv_fail=false; 88 femmodel->parameters->FindParam(&dtslices,HydrologyStepsPerStepEnum); 89 femmodel->parameters->FindParam(&sliceadapt,HydrologyStepAdaptEnum); 90 91 if(dtslices>1 || sliceadapt){ 92 int step, substep, numaveragedinput, hydro_averaging, remainingslices; 93 IssmDouble global_time, subtime, yts, remainingtime, averagetime; 91 94 IssmDouble dt, subdt; 92 95 93 96 femmodel->parameters->FindParam(&global_time,TimeEnum); 97 femmodel->parameters->FindParam(&step,StepEnum); 94 98 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); 95 99 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); … … 124 128 } 125 129 femmodel->InitTransientInputx(&transientinput[0],numaveragedinput); 130 averagetime=0; 126 131 while(substep<dtslices){ //loop on hydro dts 127 132 substep+=1; 128 133 subtime+=subdt; 134 averagetime+=subtime*dt; 129 135 /*Setting substep time as global time*/ 130 136 femmodel->parameters->SetParam(subtime,TimeEnum); … … 138 144 } 139 145 /*Proceed now to heads computations*/ 140 solutionsequence_hydro_nonlinear(femmodel); 141 /*If we have a sub-timestep we store the substep inputs in a transient input here*/ 142 femmodel->StackTransientInputonBasex(&substepinput[0],&transientinput[0],subtime,numaveragedinput); 146 solutionsequence_hydro_nonlinear(femmodel, &conv_fail); 147 if(conv_fail){ 148 /*convergence failed, we want to go back to the begining of the main step and increase the number of subslices*/ 149 /*First we get teh time and step counter back to the begining of the step that did not converge*/ 150 averagetime-=subtime*subdt; 151 subtime-=subdt; 152 substep-=1; 153 /*compute the number of slice that are remaining and the time left in the timestep*/ 154 remainingslices=dtslices-substep; 155 remainingtime=global_time-subtime; 156 /*We double the number of remaining slices and compute their duration*/ 157 dtslices=dtslices-remainingslices+(2*remainingslices); 158 subdt=remainingtime/(2*remainingslices); 159 _printf0_("convergence failed for sub-step "<< substep <<" total number of slice is now "<< dtslices <<" for step "<<step<<"\n"); 160 _printf0_("next slice duration is "<< subdt/yts <<" years\n"); 161 conv_fail = false; //re-initialize the control keyword 162 if (dtslices>500){ 163 _error_(" We reached (" << dtslices << ") which exceeds the hard limit of 500"); 164 } 165 166 femmodel->parameters->SetParam(subdt,TimesteppingTimeStepEnum); 167 femmodel->parameters->SetParam(subtime,TimeEnum); 168 169 } 170 else{ 171 /*If we have a sub-timestep we store the substep inputs in a transient input here*/ 172 femmodel->StackTransientInputonBasex(&substepinput[0],&transientinput[0],subtime,numaveragedinput); 173 } 143 174 } 144 175 /*averaging the stack*/ … … 147 178 femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum); 148 179 femmodel->parameters->SetParam(global_time,TimeEnum); 180 if(save_results){ 181 femmodel->results->AddResult(new GenericExternalResult<int>(femmodel->results->Size()+1,HydrologySubstepsEnum,dtslices,step,global_time)); 182 femmodel->results->AddResult(new GenericExternalResult<double>(femmodel->results->Size()+1,HydrologySubTimeEnum,(averagetime/dt)/yts,step,global_time)); 183 } 149 184 } 150 185 else{ … … 156 191 /*Proceed now to heads computations*/ 157 192 if(VerboseSolution()) _printf0_(" computing water heads\n"); 158 solutionsequence_hydro_nonlinear(femmodel );193 solutionsequence_hydro_nonlinear(femmodel, &conv_fail); 159 194 /*If no substeps are present we want to duplicate the results for coupling purposes*/ 160 195 InputDuplicatex(femmodel,SedimentHeadSubstepEnum,SedimentHeadEnum); -
issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
r27102 r27177 10 10 #include "../solutionsequences/solutionsequences.h" 11 11 12 void solutionsequence_hydro_nonlinear(FemModel* femmodel ){12 void solutionsequence_hydro_nonlinear(FemModel* femmodel, bool* pconv_fail){ 13 13 /*solution : */ 14 14 Vector<IssmDouble>* ug_sed=NULL; … … 16 16 Vector<IssmDouble>* uf_sed_sub_iter=NULL; 17 17 Vector<IssmDouble>* ug_sed_main_iter=NULL; 18 Vector<IssmDouble>* ug_sed_init=NULL; 18 19 19 20 Vector<IssmDouble>* ug_epl=NULL; 20 21 Vector<IssmDouble>* uf_epl=NULL; 21 22 Vector<IssmDouble>* uf_epl_sub_iter=NULL; 22 Vector<IssmDouble>* ug_epl_sub_iter=NULL;23 23 Vector<IssmDouble>* ug_epl_main_iter=NULL; 24 Vector<IssmDouble>* ug_epl_init=NULL; 24 25 25 26 Vector<IssmDouble>* ys=NULL; 26 27 Vector<IssmDouble>* dug=NULL; 27 Vector<IssmDouble>* duf=NULL;28 28 29 29 Matrix<IssmDouble>* Kff=NULL; … … 38 38 bool sedconverged,eplconverged,hydroconverged; 39 39 bool isefficientlayer; 40 bool sliceadapt; 40 41 int constraints_converged; 41 42 int num_unstable_constraints; 42 43 int sedcount,eplcount,hydrocount; 43 44 int hydro_maxiter; 45 int epl_fsize,epl_sub_fsize,epl_main_fsize; 44 46 IssmDouble sediment_kmax; 45 IssmDouble eps_ hyd;47 IssmDouble eps_res,eps_rel,eps_abs; 46 48 IssmDouble ndu_sed,nu_sed; 47 49 IssmDouble ndu_epl,nu_epl; … … 50 52 femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); 51 53 femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); 54 femmodel->parameters->FindParam(&sliceadapt,HydrologyStepAdaptEnum); 52 55 femmodel->parameters->FindParam(&hydro_maxiter,HydrologydcMaxIterEnum); 53 femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum); 56 femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum); 57 femmodel->parameters->FindParam(&eps_rel,HydrologydcRelTolEnum); 58 femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum); 54 59 hydrocount=1; 55 60 hydroconverged=false; … … 59 64 /*{{{*//*Retrieve inputs as the initial state for the non linear iteration*/ 60 65 GetBasalSolutionFromInputsx(&ug_sed,femmodel); 66 /*Initialize the IDS element mask to exclude frozen nodes*/ 67 inefanalysis->ElementizeIdsMask(femmodel); 68 61 69 Reducevectorgtofx(&uf_sed, ug_sed, femmodel->nodes,femmodel->parameters); 70 ug_sed_init=ug_sed->Duplicate(); 71 ug_sed->Copy(ug_sed_init); 72 62 73 if(isefficientlayer) { 63 74 inefanalysis = new HydrologyDCInefficientAnalysis(); … … 65 76 femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); 66 77 GetBasalSolutionFromInputsx(&ug_epl,femmodel); 67 /*Initialize the EPL element mask*/68 78 inefanalysis->ElementizeEplMask(femmodel); 69 79 effanalysis->InitZigZagCounter(femmodel); 70 /*Initialize the IDS element mask*/71 femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);72 inefanalysis->ElementizeIdsMask(femmodel);80 Reducevectorgtofx(&uf_epl, ug_epl, femmodel->nodes,femmodel->parameters); 81 ug_epl_init=ug_epl->Duplicate(); 82 ug_epl->Copy(ug_epl_init); 73 83 } 74 84 /*}}}*/ 75 85 /*The real computation starts here, outermost loop is on the two layer system*/ 76 86 for(;;){ 77 78 87 sedcount=1; 79 88 eplcount=1; … … 105 114 /*{{{*//*Core of the computation*/ 106 115 if(VerboseSolution()) _printf0_("Building Sediment Matrix...\n"); 116 107 117 femmodel->profiler->Start(SEDMatrix); 108 118 SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax,femmodel); … … 114 124 femmodel->profiler->Start(SOLVER); 115 125 Solverx(&uf_sed,Kff,pf,uf_sed_sub_iter,df,femmodel->parameters); 126 delete df; 116 127 femmodel->profiler->Stop(SOLVER); 117 118 delete Kff; delete pf; delete df;119 128 delete ug_sed; 120 129 femmodel->profiler->Start(SEDUpdate); … … 125 134 /*}}}*/ 126 135 if (!sedconverged){ 136 /*First check that all the penalizations are applied*/ 127 137 if(VerboseConvergence()) _printf0_(" # Sediment unstable constraints = " << num_unstable_constraints << "\n"); 128 if(num_unstable_constraints==0) sedconverged = true; 138 if(num_unstable_constraints==0) { 139 sedconverged = true; 140 } 141 else{//clean up 142 delete Kff; 143 delete pf; 144 } 129 145 if (sedcount>=hydro_maxiter){ 146 delete ug_sed;delete uf_sed;delete inefanalysis; delete ug_sed_main_iter; 147 if(isefficientlayer)delete ug_epl;delete uf_epl;delete effanalysis; delete ug_epl_main_iter; 130 148 _error_(" maximum number of Sediment iterations (" << hydro_maxiter << ") exceeded"); 149 131 150 } 132 151 } … … 135 154 if(sedconverged)break; 136 155 } 137 138 156 /*}}}*//*End of the sediment penalization loop*/ 139 157 sedconverged=false; 140 141 158 /*Checking convergence on the value of the sediment head*/ 142 duf=uf_sed_sub_iter->Duplicate();_assert_(duf); 143 uf_sed_sub_iter->Copy(duf); 144 duf->AYPX(uf_sed,-1.0); 145 ndu_sed=duf->Norm(NORM_TWO); 146 delete duf; 147 nu_sed=uf_sed_sub_iter->Norm(NORM_TWO); 148 if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("convergence criterion is NaN!"); 149 if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the layer is empty*/ 150 if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << "%, aiming lower than " << eps_hyd*10*100 << " %\n"); 151 if((ndu_sed/nu_sed)<eps_hyd*10.){ 152 if(VerboseConvergence()) _printf0_(" # Inner sediment convergence achieve \n"); 153 sedconverged=true; 154 } 155 delete uf_sed_sub_iter; 159 convergence(&sedconverged,Kff,pf,uf_sed,uf_sed_sub_iter,eps_res,eps_rel,eps_abs); 160 delete Kff; delete pf;delete uf_sed_sub_iter; 156 161 if(sedconverged){ 157 162 femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum); … … 176 181 for(;;){ 177 182 eplconverged=false; 178 ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter);179 ug_epl->Copy(ug_epl_sub_iter);180 183 /*{{{*//*Retrieve the EPL head slopes and compute EPL Thickness*/ 181 184 if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); … … 195 198 femmodel->profiler->Stop(EPLMasking); 196 199 if(VerboseSolution()) _printf0_("Building EPL Matrix...\n"); 200 uf_epl_sub_iter=uf_epl->Duplicate();_assert_(uf_epl_sub_iter); 201 uf_epl->Copy(uf_epl_sub_iter); 202 uf_epl->GetSize(&epl_sub_fsize); 197 203 femmodel->profiler->Start(EPLMatrices); 198 204 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); 199 205 CreateNodalConstraintsx(&ys,femmodel->nodes); 200 Reduceloadx(pf,Kfs,ys); delete Kfs;201 delete uf_epl;206 Reduceloadx(pf,Kfs,ys); 207 delete Kfs;delete uf_epl; 202 208 femmodel->profiler->Stop(EPLMatrices); 203 209 femmodel->profiler->Start(SOLVER); 204 210 Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters); 205 211 femmodel->profiler->Stop(SOLVER); 206 207 delete Kff; delete pf; delete df; 208 delete uf_epl_sub_iter; 209 uf_epl_sub_iter=uf_epl->Duplicate();_assert_(uf_epl_sub_iter); 210 uf_epl->Copy(uf_epl_sub_iter); 211 delete ug_epl; 212 delete df;delete ug_epl; 212 213 femmodel->profiler->Start(EPLUpdate); 213 214 Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys; … … 215 216 ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); 216 217 femmodel->profiler->Stop(EPLUpdate); 217 218 dug=ug_epl_sub_iter->Duplicate();_assert_(dug); 219 ug_epl_sub_iter->Copy(dug); 220 dug->AYPX(ug_epl,-1.0); 221 ndu_epl=dug->Norm(NORM_TWO); 222 delete dug; 223 nu_epl=ug_epl_sub_iter->Norm(NORM_TWO); 224 if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("convergence criterion is NaN!"); 225 if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/ 226 if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner EPL Convergence criterion:" << ndu_epl/nu_epl*100 << "%, aiming lower than " << eps_hyd*10*100 << " %\n"); 227 if((ndu_epl/nu_epl)<eps_hyd*10.) eplconverged=true; 228 if (eplcount>=hydro_maxiter){ 218 uf_epl->GetSize(&epl_fsize); 219 if(epl_fsize-epl_sub_fsize==0){ 220 convergence(&eplconverged,Kff,pf,uf_epl,uf_epl_sub_iter,eps_res,eps_rel,eps_abs); 221 delete Kff; delete pf; 222 /* if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /\*Hacking the case where the EPL is used but empty*\/ */ 223 } 224 else{ 225 delete Kff; delete pf; 226 } 227 228 if (eplcount>=hydro_maxiter*9/10 && sliceadapt && !eplconverged) { 229 if(VerboseSolution()) _printf0_("epl did not converged after "<<eplconverged<<" iteration, we refine the steping\n"); 230 *pconv_fail = true; 231 InputUpdateFromSolutionx(femmodel,ug_epl_init); 232 delete ug_epl_init; 233 femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); 234 InputUpdateFromSolutionx(femmodel,ug_sed_init); 235 delete ug_sed_init; 236 break; 237 } 238 else if (eplcount>=hydro_maxiter){ 239 delete ug_sed;delete uf_sed;delete inefanalysis;delete ug_sed_main_iter; 240 delete ug_epl;delete uf_epl;delete effanalysis;delete ug_epl_main_iter; 229 241 _error_(" maximum number of EPL iterations (" << hydro_maxiter << ") exceeded"); 242 230 243 } 231 244 eplcount++; 232 233 delete ug_epl_sub_iter; 245 delete uf_epl_sub_iter; 234 246 if(eplconverged){ 235 247 if(VerboseSolution()) _printf0_("eplconverged...\n"); 236 InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum);237 InputUpdateFromSolutionx(femmodel,ug_epl);238 248 effanalysis->ResetCounter(femmodel); 239 249 break; … … 242 252 } 243 253 femmodel->profiler->Stop(EPLLOOP); 254 if(*pconv_fail){ 255 break; 256 } 244 257 /*}}}*//*End of the global EPL loop*/ 245 258 /*{{{*//*Now dealing with the convergence of the whole system*/ … … 264 277 if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("EPL convergence criterion is NaN!"); 265 278 if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/ 266 if (!xIsNan<IssmDouble>(eps_ hyd)){267 if ((ndu_epl/nu_epl)<eps_ hyd && (ndu_sed/nu_sed)<(eps_hyd)){279 if (!xIsNan<IssmDouble>(eps_rel)){ 280 if ((ndu_epl/nu_epl)<eps_rel && (ndu_sed/nu_sed)<(eps_rel)){ 268 281 if (VerboseConvergence()) _printf0_(setw(50) << left << " Converged after, " << hydrocount << " iterations \n"); 269 282 hydroconverged=true; 270 283 } 271 284 else{ 272 if(VerboseConvergence()) _printf0_(setw(50) << left << " Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << "%, aiming lower than " << eps_hyd*100 << " %\n"); 273 if(VerboseConvergence()) _printf0_(setw(50) << left << " EPL Convergence criterion:" << ndu_epl/nu_epl*100 << "%, aiming lower than " << eps_hyd*100 << " %\n"); 285 if(VerboseConvergence()) _printf0_(setw(50) << left << " for iteration:" << hydrocount << " \n"); 286 if(VerboseConvergence()) _printf0_(setw(50) << left << " Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << "%, aiming lower than " << eps_rel*100 << " %\n"); 287 if(VerboseConvergence()) _printf0_(setw(50) << left << " EPL Convergence criterion:" << ndu_epl/nu_epl*100 << "%, aiming lower than " << eps_rel*100 << " %\n"); 274 288 hydroconverged=false; 275 289 } 276 290 } 277 291 else _printf0_(setw(50) << left << " Convergence criterion:" << ndu_sed/nu_sed*100 << " %\n"); 292 293 if (hydrocount>=hydro_maxiter*9/10 && sliceadapt && !hydroconverged) { 294 if(VerboseSolution()) _printf0_("hydrology main loop did not converged after "<<hydrocount<<" iteration, we refine the steping\n"); 295 *pconv_fail = true; 296 InputUpdateFromSolutionx(femmodel,ug_epl_init); 297 delete ug_epl_init; 298 femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); 299 InputUpdateFromSolutionx(femmodel,ug_sed_init); 300 delete ug_sed_init; 301 break; 302 } 278 303 if (hydrocount>=hydro_maxiter){ 279 304 _error_(" maximum number for hydrological global iterations (" << hydro_maxiter << ") exceeded"); 305 delete ug_sed;delete uf_sed;delete effanalysis; 306 delete ug_epl; delete uf_epl; delete inefanalysis; 280 307 } 281 308 } … … 284 311 } 285 312 /*}}}*/ 286 if(isefficientlayer)InputUpdateFromSolutionx(femmodel,ug_epl); 287 femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); 288 InputUpdateFromSolutionx(femmodel,ug_sed); 313 /*To deal with adaptative stepping we only save results if we are actually converged*/ 314 if(hydroconverged){ 315 if(isefficientlayer){ 316 delete ug_epl_init; 317 } 318 delete ug_sed_init; 319 } 289 320 /*Free resources: */ 290 delete ug_epl; 291 delete ug_sed; 292 delete uf_sed; 293 delete uf_epl; 294 delete uf_epl_sub_iter; 295 delete inefanalysis; 296 delete effanalysis; 321 delete ug_epl;delete ug_sed; 322 delete uf_sed; delete uf_epl; 323 delete inefanalysis; delete effanalysis; 297 324 } -
issm/trunk-jpl/src/m/classes/hydrologydc.m
r25024 r27177 14 14 max_iter = 0; 15 15 steps_per_step = 0; 16 step_adapt = 0; 16 17 averaging = 0; 17 18 sedimentlimit_flag = 0; … … 66 67 list=[list,{'EplHead','HydrologydcMaskEplactiveNode','HydrologydcMaskEplactiveElt','EplHeadSlopeX','EplHeadSlopeY','HydrologydcEplThickness'}]; 67 68 end 68 if self.steps_per_step>1 ,69 if self.steps_per_step>1 | self.step_adapt, 69 70 list = [list,'EffectivePressureSubstep','SedimentHeadSubstep']; 70 71 if self.isefficientlayer, … … 90 91 self.max_iter = 100; 91 92 self.steps_per_step = 1; 93 self.step_adapt = 0; 92 94 self.averaging = 0; 93 95 self.sedimentlimit_flag = 0; … … 126 128 md = checkfield(md,'fieldname','hydrology.max_iter','>',0,'numel',1); 127 129 md = checkfield(md,'fieldname','hydrology.steps_per_step','>',0,'numel',1); 130 md = checkfield(md,'fieldname','hydrology.step_adapt','numel',1,'values',[0 1]); 128 131 md = checkfield(md,'fieldname','hydrology.averaging','numel',[1],'values',[0 1 2]); 129 132 md = checkfield(md,'fieldname','hydrology.sedimentlimit_flag','numel',[1],'values',[0 1 2 3]); … … 172 175 fielddisplay(self,'max_iter','maximum number of nonlinear iteration'); 173 176 fielddisplay(self,'steps_per_step','number of hydrology steps per timestep'); 177 fielddisplay(self,'step_adapt', 'adaptative sub stepping [1: true 0: false] default is 0'); 174 178 fielddisplay(self, 'averaging', 'averaging methods from short to long steps'); 175 179 disp(sprintf('%55s 0: Arithmetic (default)')); … … 226 230 WriteData(fid,prefix,'object',self,'fieldname','max_iter','format','Integer'); 227 231 WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer'); 232 WriteData(fid,prefix,'object',self,'fieldname','step_adapt','format','Boolean'); 228 233 WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer'); 229 234 WriteData(fid,prefix,'object',self,'fieldname','sedimentlimit_flag','format','Integer'); -
issm/trunk-jpl/src/m/classes/hydrologydc.py
r26059 r27177 21 21 self.max_iter = 0 22 22 self.steps_per_step = 0 23 self.step_adapt = 0 23 24 self.averaging = 0 24 25 self.sedimentlimit_flag = 0 … … 53 54 def __repr__(self): # {{{ 54 55 # TODO: 55 # - Convert all formatting to calls to <string>.format (see any 56 # - Convert all formatting to calls to <string>.format (see any 56 57 # already converted <class>.__repr__ method for examples) 57 58 # … … 65 66 string = "%s\n%s" % (string, fielddisplay(self, 'max_iter', 'maximum number of nonlinear iteration')) 66 67 string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of hydrology steps per time step')) 68 string = "%s\n%s" % (string, fielddisplay(self, 'step_adapt', 'adaptative sub stepping [1: true 0: false] default is 0')) 67 69 string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps')) 68 70 string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)') … … 134 136 self.max_iter = 100 135 137 self.steps_per_step = 1 138 self.step_adapt = 0 136 139 self.averaging = 0 137 140 self.sedimentlimit_flag = 0 … … 162 165 if self.isefficientlayer == 1: 163 166 list.extend(['EplHead', 'HydrologydcMaskEplactiveNode', 'HydrologydcMaskEplactiveElt', 'EplHeadSlopeX', 'EplHeadSlopeY', 'HydrologydcEplThickness']) 164 if self.steps_per_step > 1 :167 if self.steps_per_step > 1 or self.step_adapt: 165 168 list.extend(['EffectivePressureSubstep', 'SedimentHeadSubstep']) 166 169 if self.isefficientlayer == 1: … … 190 193 md = checkfield(md, 'fieldname', 'hydrology.max_iter', '>', 0., 'numel', [1]) 191 194 md = checkfield(md, 'fieldname', 'hydrology.steps_per_step', '>=', 1, 'numel', [1]) 195 md = checkfield(md, 'fieldname', 'hydrology.step_adapt', 'numel', [1], 'values', [0, 1]) 192 196 md = checkfield(md, 'fieldname', 'hydrology.averaging', 'numel', [1], 'values', [0, 1, 2]) 193 197 md = checkfield(md, 'fieldname', 'hydrology.sedimentlimit_flag', 'numel', [1], 'values', [0, 1, 2, 3]) … … 233 237 WriteData(fid, prefix, 'object', self, 'fieldname', 'max_iter', 'format', 'Integer') 234 238 WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer') 239 WriteData(fid, prefix, 'object', self, 'fieldname', 'step_adapt', 'format', 'Boolean') 235 240 WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer') 236 241 WriteData(fid, prefix, 'object', self, 'fieldname', 'sedimentlimit_flag', 'format', 'Integer')
Note:
See TracChangeset
for help on using the changeset viewer.
