Ice Sheet System Model  4.18
Code documentation
solutionsequence_hydro_nonlinear.cpp
Go to the documentation of this file.
1 
5 #include "../toolkits/toolkits.h"
6 #include "../classes/classes.h"
7 #include "../shared/shared.h"
8 #include "../modules/modules.h"
9 /*FIXME, dirty hack to get the solutionsequence linear needed to compute the slopes*/
10 #include "../solutionsequences/solutionsequences.h"
11 
12 void solutionsequence_hydro_nonlinear(FemModel* femmodel){
13  /*solution : */
14  Vector<IssmDouble>* ug_sed=NULL;
15  Vector<IssmDouble>* uf_sed=NULL;
16  Vector<IssmDouble>* uf_sed_sub_iter=NULL;
17  Vector<IssmDouble>* ug_sed_main_iter=NULL;
18 
19  Vector<IssmDouble>* ug_epl=NULL;
20  Vector<IssmDouble>* uf_epl=NULL;
21  Vector<IssmDouble>* uf_epl_sub_iter=NULL;
22  Vector<IssmDouble>* ug_epl_sub_iter=NULL;
23  Vector<IssmDouble>* ug_epl_main_iter=NULL;
24 
25  Vector<IssmDouble>* ys=NULL;
26  Vector<IssmDouble>* dug=NULL;
27  Vector<IssmDouble>* duf=NULL;
28 
29  Matrix<IssmDouble>* Kff=NULL;
30  Matrix<IssmDouble>* Kfs=NULL;
31 
32  Vector<IssmDouble>* pf=NULL;
33  Vector<IssmDouble>* df=NULL;
34 
35  HydrologyDCInefficientAnalysis* inefanalysis = NULL;
36  HydrologyDCEfficientAnalysis* effanalysis = NULL;
37 
38  bool sedconverged,eplconverged,hydroconverged;
39  bool isefficientlayer;
40  int constraints_converged;
41  int num_unstable_constraints;
42  int sedcount,eplcount,hydrocount;
43  int hydro_maxiter;
44  IssmDouble sediment_kmax;
45  IssmDouble eps_hyd;
46  IssmDouble ndu_sed,nu_sed;
47  IssmDouble ndu_epl,nu_epl;
48  IssmDouble ThickCount,L2Count;
49  /*Recover parameters: */
50  femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
51  femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
52  femmodel->parameters->FindParam(&hydro_maxiter,HydrologydcMaxIterEnum);
53  femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum);
54  hydrocount=1;
55  hydroconverged=false;
56  /*We don't need the outer loop if only one layer is used*/
57  if(!isefficientlayer) hydroconverged=true;
58 
59  /*{{{*//*Retrieve inputs as the initial state for the non linear iteration*/
60  GetSolutionFromInputsx(&ug_sed,femmodel);
61  Reducevectorgtofx(&uf_sed, ug_sed, femmodel->nodes,femmodel->parameters);
62  if(isefficientlayer) {
63  inefanalysis = new HydrologyDCInefficientAnalysis();
64  effanalysis = new HydrologyDCEfficientAnalysis();
65  femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
66  GetSolutionFromInputsx(&ug_epl,femmodel);
67  /*Initialize the EPL element mask*/
68  inefanalysis->ElementizeEplMask(femmodel);
69  effanalysis->InitZigZagCounter(femmodel);
70  /*Initialize the IDS element mask*/
71  femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
72  inefanalysis->ElementizeIdsMask(femmodel);
73  }
74  /*}}}*/
75  /*The real computation starts here, outermost loop is on the two layer system*/
76  for(;;){
77 
78  sedcount=1;
79  eplcount=1;
80 
81  /*If there is two layers we need an outer loop value to compute convergence*/
82  if(isefficientlayer){
83  ug_sed_main_iter=ug_sed->Duplicate();
84  ug_sed->Copy(ug_sed_main_iter);
85  ug_epl_main_iter=ug_epl->Duplicate();
86  ug_epl->Copy(ug_epl_main_iter);
87  }
88  /*Loop on sediment layer to deal with transfer and head value*/
89  femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
90  ResetZigzagCounterx(femmodel);
91  InputUpdateFromConstantx(femmodel,false,ConvergedEnum);
92  femmodel->UpdateConstraintsx();
93 
94  /*Reset constraint on the ZigZag Lock*/
95  ResetConstraintsx(femmodel);
96 
97  /*{{{*//*Treating the sediment*/
98  for(;;){
99  sedconverged=false;
100  uf_sed_sub_iter=uf_sed->Duplicate();_assert_(uf_sed_sub_iter);
101  uf_sed->Copy(uf_sed_sub_iter);
102  /*{{{*//*Loop on the sediment layer to deal with the penalization*/
103  for(;;){
104  /*{{{*//*Core of the computation*/
105  if(VerboseSolution()) _printf0_("Building Sediment Matrix...\n");
106  SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax,femmodel);
107  CreateNodalConstraintsx(&ys,femmodel->nodes);
108  Reduceloadx(pf,Kfs,ys); delete Kfs;
109  delete uf_sed;
110 
111  femmodel->profiler->Start(SOLVER);
112  Solverx(&uf_sed,Kff,pf,uf_sed_sub_iter,df,femmodel->parameters);
113  femmodel->profiler->Stop(SOLVER);
114 
115  delete Kff; delete pf; delete df;
116  delete ug_sed;
117  Mergesolutionfromftogx(&ug_sed,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys;
118  InputUpdateFromSolutionx(femmodel,ug_sed);
119  ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
120  /*}}}*/
121  if (!sedconverged){
122  if(VerboseConvergence()) _printf0_(" # Sediment unstable constraints = " << num_unstable_constraints << "\n");
123  if(num_unstable_constraints==0) sedconverged = true;
124  if (sedcount>=hydro_maxiter){
125  _error_(" maximum number of Sediment iterations (" << hydro_maxiter << ") exceeded");
126  }
127  }
128  /*Add an iteration and get out of the loop if the penalisation is converged*/
129  sedcount++;
130  if(sedconverged)break;
131  }
132 
133  /*}}}*//*End of the sediment penalization loop*/
134  sedconverged=false;
135 
136  /*Checking convergence on the value of the sediment head*/
137  duf=uf_sed_sub_iter->Duplicate();_assert_(duf);
138  uf_sed_sub_iter->Copy(duf);
139  duf->AYPX(uf_sed,-1.0);
140  ndu_sed=duf->Norm(NORM_TWO);
141  delete duf;
142  nu_sed=uf_sed_sub_iter->Norm(NORM_TWO);
143  if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("convergence criterion is NaN!");
144  if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the layer is empty*/
145  if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << "%, aiming lower than " << eps_hyd*10*100 << " %\n");
146  if((ndu_sed/nu_sed)<eps_hyd*10.){
147  if(VerboseConvergence()) _printf0_(" # Inner sediment convergence achieve \n");
148  sedconverged=true;
149  }
150  delete uf_sed_sub_iter;
151  if(sedconverged){
152  femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum);
153  InputUpdateFromConstantx(femmodel,sedconverged,ConvergedEnum);
154  InputUpdateFromSolutionx(femmodel,ug_sed);
155  InputUpdateFromConstantx(femmodel,sediment_kmax,HydrologySedimentKmaxEnum);
156  break;
157  }
158  }
159  /*}}}*//*End of the global sediment loop*/
160  /*{{{*//*Now dealing with the EPL in the same way*/
161  if(isefficientlayer){
162  femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
163  /*updating mask*/
164  if(VerboseSolution()) _printf0_("==updating mask...\n");
165  femmodel->HydrologyEPLupdateDomainx(&ThickCount);
166  ResetZigzagCounterx(femmodel);
167  InputUpdateFromConstantx(femmodel,false,ConvergedEnum);
168 
169  for(;;){
170  eplconverged=false;
171  ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter);
172  ug_epl->Copy(ug_epl_sub_iter);
173  /*{{{*//*Retrieve the EPL head slopes and compute EPL Thickness*/
174  if(VerboseSolution()) _printf0_("computing EPL Head slope...\n");
175  femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum);
176  femmodel->UpdateConstraintsL2ProjectionEPLx(&L2Count);
177  femmodel->parameters->SetParam(EplHeadSlopeXEnum,InputToL2ProjectEnum);
178  solutionsequence_linear(femmodel);
179  femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum);
180  solutionsequence_linear(femmodel);
181 
182  femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
183  effanalysis->ComputeEPLThickness(femmodel);
184  //updating mask after the computation of the epl thickness (Allow to close too thin EPL)
185  femmodel->HydrologyEPLupdateDomainx(&ThickCount);
186  /*}}}*/
187 
188  if(VerboseSolution()) _printf0_("Building EPL Matrix...\n");
189  SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
190  CreateNodalConstraintsx(&ys,femmodel->nodes);
191  Reduceloadx(pf,Kfs,ys); delete Kfs;
192  delete uf_epl;
193 
194  femmodel->profiler->Start(SOLVER);
195  Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters);
196  femmodel->profiler->Stop(SOLVER);
197 
198  delete Kff; delete pf; delete df;
199  delete uf_epl_sub_iter;
200  uf_epl_sub_iter=uf_epl->Duplicate();_assert_(uf_epl_sub_iter);
201  uf_epl->Copy(uf_epl_sub_iter);
202  delete ug_epl;
203  Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys;
204  InputUpdateFromSolutionx(femmodel,ug_epl);
205  ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
206 
207  dug=ug_epl_sub_iter->Duplicate();_assert_(dug);
208  ug_epl_sub_iter->Copy(dug);
209  dug->AYPX(ug_epl,-1.0);
210  ndu_epl=dug->Norm(NORM_TWO);
211  delete dug;
212  nu_epl=ug_epl_sub_iter->Norm(NORM_TWO);
213  if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("convergence criterion is NaN!");
214  if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/
215  if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner EPL Convergence criterion:" << ndu_epl/nu_epl*100 << "%, aiming lower than " << eps_hyd*10*100 << " %\n");
216  if((ndu_epl/nu_epl)<eps_hyd*10.) eplconverged=true;
217  if (eplcount>=hydro_maxiter){
218  _error_(" maximum number of EPL iterations (" << hydro_maxiter << ") exceeded");
219  }
220  eplcount++;
221 
222  delete ug_epl_sub_iter;
223  if(eplconverged){
224  if(VerboseSolution()) _printf0_("eplconverged...\n");
225  InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum);
226  InputUpdateFromSolutionx(femmodel,ug_epl);
227  effanalysis->ResetCounter(femmodel);
228  break;
229  }
230  }
231  }
232  /*}}}*//*End of the global EPL loop*/
233  /*{{{*//*Now dealing with the convergence of the whole system*/
234  if(!hydroconverged){
235  //compute norm(du)/norm(u)
236  dug=ug_sed_main_iter->Duplicate(); _assert_(dug);
237  ug_sed_main_iter->Copy(dug);
238  dug->AYPX(ug_sed,-1.0);
239  ndu_sed=dug->Norm(NORM_TWO);
240  delete dug;
241  nu_sed=ug_sed_main_iter->Norm(NORM_TWO);
242  delete ug_sed_main_iter;
243  if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("Sed convergence criterion is NaN!");
244  if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the Sediment is used but empty*/
245  dug=ug_epl_main_iter->Duplicate();_assert_(dug);
246  ug_epl_main_iter->Copy(dug);
247  dug->AYPX(ug_epl,-1.0);
248  ndu_epl=dug->Norm(NORM_TWO);
249  delete dug;
250  nu_epl=ug_epl_main_iter->Norm(NORM_TWO);
251  delete ug_epl_main_iter;
252  if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("EPL convergence criterion is NaN!");
253  if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/
254  if (!xIsNan<IssmDouble>(eps_hyd)){
255  if ((ndu_epl/nu_epl)<eps_hyd && (ndu_sed/nu_sed)<(eps_hyd)){
256  if (VerboseConvergence()) _printf0_(setw(50) << left << " Converged after, " << hydrocount << " iterations \n");
257  hydroconverged=true;
258  }
259  else{
260  if(VerboseConvergence()) _printf0_(setw(50) << left << " Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << "%, aiming lower than " << eps_hyd*100 << " %\n");
261  if(VerboseConvergence()) _printf0_(setw(50) << left << " EPL Convergence criterion:" << ndu_epl/nu_epl*100 << "%, aiming lower than " << eps_hyd*100 << " %\n");
262  hydroconverged=false;
263  }
264  }
265  else _printf0_(setw(50) << left << " Convergence criterion:" << ndu_sed/nu_sed*100 << " %\n");
266  if (hydrocount>=hydro_maxiter){
267  _error_(" maximum number for hydrological global iterations (" << hydro_maxiter << ") exceeded");
268  }
269  }
270  hydrocount++;
271  if(hydroconverged)break;
272  }
273  /*}}}*/
274  if(isefficientlayer)InputUpdateFromSolutionx(femmodel,ug_epl);
275  femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
276  InputUpdateFromSolutionx(femmodel,ug_sed);
277  /*Free ressources: */
278  delete ug_epl;
279  delete ug_sed;
280  delete uf_sed;
281  delete uf_epl;
282  delete uf_epl_sub_iter;
283  delete inefanalysis;
284  delete effanalysis;
285 }
Matrix< IssmDouble >
SOLVER
#define SOLVER
Definition: Profiler.h:16
Mergesolutionfromftogx
void Mergesolutionfromftogx(Vector< IssmDouble > **pug, Vector< IssmDouble > *uf, Vector< IssmDouble > *ys, Nodes *nodes, Parameters *parameters, bool flag_ys0)
Definition: Mergesolutionfromftogx.cpp:8
EplHeadSlopeXEnum
@ EplHeadSlopeXEnum
Definition: EnumDefinitions.h:555
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
ConstraintsStatex
void ConstraintsStatex(int *pconverged, int *pnum_unstable_constraints, FemModel *femmodel)
Definition: ConstraintsStatex.cpp:10
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
VerboseConvergence
bool VerboseConvergence(void)
Definition: Verbosity.cpp:26
ResetZigzagCounterx
void ResetZigzagCounterx(FemModel *femmodel)
Definition: ResetConstraintsx.cpp:39
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
ConvergedEnum
@ ConvergedEnum
Definition: EnumDefinitions.h:514
HydrologydcIsefficientlayerEnum
@ HydrologydcIsefficientlayerEnum
Definition: EnumDefinitions.h:185
Solverx
void Solverx(Vector< IssmDouble > **puf, Matrix< IssmDouble > *Kff, Vector< IssmDouble > *pf, Vector< IssmDouble > *uf0, Vector< IssmDouble > *df, Parameters *parameters)
Definition: Solverx.cpp:15
Vector::Duplicate
Vector< doubletype > * Duplicate(void)
Definition: Vector.h:230
HydrologydcRelTolEnum
@ HydrologydcRelTolEnum
Definition: EnumDefinitions.h:190
solutionsequence_linear
void solutionsequence_linear(FemModel *femmodel)
Definition: solutionsequence_linear.cpp:10
NORM_TWO
@ NORM_TWO
Definition: toolkitsenums.h:15
HydrologyDCEfficientAnalysis
Definition: HydrologyDCEfficientAnalysis.h:12
HydrologyDCInefficientAnalysisEnum
@ HydrologyDCInefficientAnalysisEnum
Definition: EnumDefinitions.h:1099
Vector::Norm
doubletype Norm(NormMode norm_type)
Definition: Vector.h:342
HydrologyDCEfficientAnalysis::InitZigZagCounter
void InitZigZagCounter(FemModel *femmodel)
Definition: HydrologyDCEfficientAnalysis.cpp:136
Parameters::SetParam
void SetParam(bool boolean, int enum_type)
Definition: Parameters.cpp:441
HydrologySedimentKmaxEnum
@ HydrologySedimentKmaxEnum
Definition: EnumDefinitions.h:174
FemModel::HydrologyEPLupdateDomainx
void HydrologyEPLupdateDomainx(IssmDouble *pEplcount)
Definition: FemModel.cpp:4965
FemModel::UpdateConstraintsx
void UpdateConstraintsx(void)
Definition: FemModel.cpp:3027
FemModel::nodes
Nodes * nodes
Definition: FemModel.h:56
HydrologyDCInefficientAnalysis::ElementizeIdsMask
void ElementizeIdsMask(FemModel *femmodel)
Definition: HydrologyDCInefficientAnalysis.cpp:816
HydrologyDCInefficientAnalysis
Definition: HydrologyDCInefficientAnalysis.h:12
HydrologydcMaxIterEnum
@ HydrologydcMaxIterEnum
Definition: EnumDefinitions.h:187
HydrologyDCEfficientAnalysis::ComputeEPLThickness
void ComputeEPLThickness(FemModel *femmodel)
Definition: HydrologyDCEfficientAnalysis.cpp:615
InputUpdateFromConstantx
void InputUpdateFromConstantx(FemModel *femmodel, bool constant, int name)
Definition: InputUpdateFromConstantx.cpp:10
FemModel::UpdateConstraintsL2ProjectionEPLx
void UpdateConstraintsL2ProjectionEPLx(IssmDouble *pL2count)
Definition: FemModel.cpp:5160
Profiler::Stop
void Stop(int tagenum, bool dontmpisync=true)
Definition: Profiler.cpp:179
HydrologyDCEfficientAnalysisEnum
@ HydrologyDCEfficientAnalysisEnum
Definition: EnumDefinitions.h:1098
Reducevectorgtofx
void Reducevectorgtofx(Vector< IssmDouble > **puf, Vector< IssmDouble > *ug, Nodes *nodes, Parameters *parameters)
Definition: Reducevectorgtofx.cpp:8
FemModel
Definition: FemModel.h:31
Reduceloadx
void Reduceloadx(Vector< IssmDouble > *pf, Matrix< IssmDouble > *Kfs, Vector< IssmDouble > *y_s, bool flag_ys0)
Definition: Reduceloadx.cpp:14
L2ProjectionEPLAnalysisEnum
@ L2ProjectionEPLAnalysisEnum
Definition: EnumDefinitions.h:1137
Vector::AYPX
void AYPX(Vector *X, doubletype a)
Definition: Vector.h:267
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Profiler::Start
void Start(int tagenum, bool dontmpisync=true)
Definition: Profiler.cpp:139
InputUpdateFromSolutionx
void InputUpdateFromSolutionx(FemModel *femmodel, Vector< IssmDouble > *solution)
Definition: InputUpdateFromSolutionx.cpp:9
GetSolutionFromInputsx
void GetSolutionFromInputsx(Vector< IssmDouble > **psolution, FemModel *femmodel)
Definition: GetSolutionFromInputsx.cpp:9
VerboseSolution
bool VerboseSolution(void)
Definition: Verbosity.cpp:24
HydrologyDCEfficientAnalysis::ResetCounter
void ResetCounter(FemModel *femmodel)
Definition: HydrologyDCEfficientAnalysis.cpp:143
FemModel::SetCurrentConfiguration
void SetCurrentConfiguration(int configuration_type)
Definition: FemModel.cpp:634
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
ResetConstraintsx
void ResetConstraintsx(FemModel *femmodel)
Definition: ResetConstraintsx.cpp:16
InputToL2ProjectEnum
@ InputToL2ProjectEnum
Definition: EnumDefinitions.h:206
solutionsequence_hydro_nonlinear
void solutionsequence_hydro_nonlinear(FemModel *femmodel)
Definition: solutionsequence_hydro_nonlinear.cpp:12
Vector::Copy
void Copy(Vector *to)
Definition: Vector.h:319
HydrologyDCInefficientAnalysis::ElementizeEplMask
void ElementizeEplMask(FemModel *femmodel)
Definition: HydrologyDCInefficientAnalysis.cpp:756
CreateNodalConstraintsx
void CreateNodalConstraintsx(Vector< IssmDouble > **pys, Nodes *nodes)
Definition: CreateNodalConstraintsx.cpp:10
EplHeadSlopeYEnum
@ EplHeadSlopeYEnum
Definition: EnumDefinitions.h:556
SystemMatricesx
void SystemMatricesx(Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, Vector< IssmDouble > **ppf, Vector< IssmDouble > **pdf, IssmDouble *pkmax, FemModel *femmodel, bool isAllocated)
Definition: SystemMatricesx.cpp:10
FemModel::profiler
Profiler * profiler
Definition: FemModel.h:42
Vector< IssmDouble >
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:22 for Ice Sheet System Model by   doxygen 1.8.19