Ice Sheet System Model  4.18
Code documentation
Functions
solutionsequences.h File Reference
#include "../shared/Numerics/types.h"

Go to the source code of this file.

Functions

void solutionsequence_thermal_nonlinear (FemModel *femmodel)
 
void solutionsequence_hydro_nonlinear (FemModel *femmodel)
 
void solutionsequence_shakti_nonlinear (FemModel *femmodel)
 
void solutionsequence_glads_nonlinear (FemModel *femmodel)
 
void solutionsequence_nonlinear (FemModel *femmodel, bool conserve_loads)
 
void solutionsequence_newton (FemModel *femmodel)
 
void solutionsequence_fct (FemModel *femmodel)
 
void solutionsequence_FScoupling_nonlinear (FemModel *femmodel, bool conserve_loads)
 
void solutionsequence_linear (FemModel *femmodel)
 
void solutionsequence_la (FemModel *femmodel)
 
void solutionsequence_la_theta (FemModel *femmodel)
 
void solutionsequence_adjoint_linear (FemModel *femmodel)
 
void solutionsequence_schurcg (FemModel *femmodel)
 
void convergence (bool *pconverged, Matrix< IssmDouble > *K_ff, Vector< IssmDouble > *p_f, Vector< IssmDouble > *u_f, Vector< IssmDouble > *u_f_old, IssmDouble eps_res, IssmDouble eps_rel, IssmDouble eps_abs)
 

Function Documentation

◆ solutionsequence_thermal_nonlinear()

void solutionsequence_thermal_nonlinear ( FemModel *  femmodel)

Definition at line 11 of file solutionsequence_thermal_nonlinear.cpp.

11  {
12 
13  /*solution : */
14  Vector<IssmDouble>* tg=NULL;
15  Vector<IssmDouble>* tf=NULL;
16  Vector<IssmDouble>* tf_old=NULL;
17  Vector<IssmDouble>* ys=NULL;
18  IssmDouble melting_offset;
19 
20  /*intermediary: */
21  Matrix<IssmDouble>* Kff=NULL;
22  Matrix<IssmDouble>* Kfs=NULL;
23  Vector<IssmDouble>* pf=NULL;
24  Vector<IssmDouble>* df=NULL;
25 
26  bool converged;
27  bool isenthalpy, isdynamicbasalspc;
28  int constraints_converged;
29  int num_unstable_constraints;
30  int count;
31  int thermal_penalty_threshold;
32  int thermal_maxiter;
33 
34  /*parameters:*/
35  IssmDouble eps_rel;
36 
37  /*Recover parameters: */
38  femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
39  femmodel->parameters->FindParam(&thermal_maxiter,ThermalMaxiterEnum);
40 
41  converged=false;
42  InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
43 
44  if(isenthalpy){
45  femmodel->parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
46  femmodel->parameters->FindParam(&eps_rel,ThermalReltolEnum);
47  femmodel->UpdateConstraintsx();
48 
49  //Update the solution to make sure that tf and tf_old are similar (for next step in transient or steadystate)
50  GetSolutionFromInputsx(&tg,femmodel);
51  Reducevectorgtofx(&tf, tg, femmodel->nodes,femmodel->parameters);
52  InputUpdateFromSolutionx(femmodel,tg);
53  }
54  else{
55  femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum);
56  ResetZigzagCounterx(femmodel);
57  femmodel->UpdateConstraintsx();
58  }
59 
60  count=1;
61 
62  for(;;){
63  delete tf_old;tf_old=tf;
64 
65  if(isenthalpy){
66  SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
67  /*Update old solution, such that sizes of tf_old and tf are comparable*/
68  if(isdynamicbasalspc){
69  delete tf_old;
70  Reducevectorgtofx(&tf_old, tg, femmodel->nodes,femmodel->parameters);
71  }
72  }
73  else SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel);
74  delete tg;
75  CreateNodalConstraintsx(&ys,femmodel->nodes);
76  Reduceloadx(pf, Kfs, ys); delete Kfs;
77  femmodel->profiler->Start(SOLVER);
78  Solverx(&tf, Kff, pf, tf_old, df, femmodel->parameters);
79  femmodel->profiler->Stop(SOLVER);
80  Mergesolutionfromftogx(&tg, tf,ys,femmodel->nodes,femmodel->parameters); delete ys;
81  if(isenthalpy){
82  convergence(&converged,Kff,pf,tf,tf_old,0.05,eps_rel,NAN);
83  InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
84  }
85  delete Kff; delete pf; delete df;
86  InputUpdateFromSolutionx(femmodel,tg);
87  ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
88  if(VerboseConvergence()) _printf0_(" number of unstable constraints: " << num_unstable_constraints << "\n");
89 
90  if(isenthalpy){ // enthalpy method
91  IssmDouble dt;
92  femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
93 
94  count++;
95  bool max_iteration_state=false;
96  if(count>=thermal_maxiter){
97  _printf0_(" maximum number of nonlinear iterations (" << thermal_maxiter << ") exceeded\n");
98  converged=true;
99  InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
100  InputUpdateFromSolutionx(femmodel,tg);
101  max_iteration_state=true;
102  }
103  if(converged==true){
104  break;
105  }
106  else if(dt==0.){
107  EnthalpyAnalysis::ComputeBasalMeltingrate(femmodel);
108  EnthalpyAnalysis::UpdateBasalConstraints(femmodel);
109  }
110  }
111  else{ // dry ice method
112  if(!converged){
113  if(num_unstable_constraints<=thermal_penalty_threshold) converged=true;
114  if(count>=thermal_maxiter){
115  converged=true;
116  _printf0_(" maximum number of iterations (" << thermal_maxiter << ") exceeded\n");
117  }
118  }
119  count++;
120  InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
121  if(converged)break;
122  }
123  }
124 
125  if(isenthalpy){
126  if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n");
127  }
128  else{
129  InputUpdateFromSolutionx(femmodel,tg);
130  femmodel->parameters->SetParam(melting_offset,MeltingOffsetEnum);
131  }
132 
133  /*Free ressources: */
134  delete tg;
135  delete tf;
136  delete tf_old;
137 }

◆ solutionsequence_hydro_nonlinear()

void solutionsequence_hydro_nonlinear ( FemModel *  femmodel)

Definition at line 12 of file solutionsequence_hydro_nonlinear.cpp.

12  {
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 }

◆ solutionsequence_shakti_nonlinear()

void solutionsequence_shakti_nonlinear ( FemModel *  femmodel)

Definition at line 11 of file solutionsequence_shakti_nonlinear.cpp.

11  {
12 
13  /*intermediary: */
14  Matrix<IssmDouble>* Kff = NULL;
15  Matrix<IssmDouble>* Kfs = NULL;
16  Vector<IssmDouble>* ug = NULL;
17  Vector<IssmDouble>* uf = NULL;
18  Vector<IssmDouble>* old_uf = NULL;
19  Vector<IssmDouble>* pf = NULL;
20  Vector<IssmDouble>* df = NULL;
21  Vector<IssmDouble>* ys = NULL;
22 
23  /*parameters:*/
24  int max_nonlinear_iterations;
25  IssmDouble eps_res,eps_rel,eps_abs;
26 
27  /*Recover parameters: */
28  femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
29  femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum);
30  femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
31  femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
32  femmodel->UpdateConstraintsx();
33 
34  int count=0;
35  bool converged=false;
36 
37  /*Start non-linear iteration using input velocity: */
38  GetSolutionFromInputsx(&ug,femmodel);
39  Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
40 
41  for(;;){
42  /*save pointer to old solution*/
43  delete old_uf;old_uf=uf;
44  delete ug;
45 
46  SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
47  CreateNodalConstraintsx(&ys,femmodel->nodes);
48  Reduceloadx(pf, Kfs, ys); delete Kfs;
49  femmodel->profiler->Start(SOLVER);
50  Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters);
51  femmodel->profiler->Stop(SOLVER);
52  Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
53 
54  convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs); delete Kff; delete pf; delete df;
55  InputUpdateFromSolutionx(femmodel,ug);
56 
57  /*Increase count: */
58  count++;
59  if(count>=max_nonlinear_iterations){
60  _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");
61  converged = true;
62  }
63  if(converged) break;
64  }
65 
66  if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count << "\n");
67 
68  /*clean-up*/
69  delete uf;
70  delete ug;
71  delete old_uf;
72 }

◆ solutionsequence_glads_nonlinear()

void solutionsequence_glads_nonlinear ( FemModel *  femmodel)

Definition at line 11 of file solutionsequence_glads_nonlinear.cpp.

11  {
12 
13  /*intermediary: */
14  Matrix<IssmDouble>* Kff = NULL;
15  Matrix<IssmDouble>* Kfs = NULL;
16  Vector<IssmDouble>* ug = NULL;
17  Vector<IssmDouble>* uf = NULL;
18  Vector<IssmDouble>* old_uf = NULL;
19  Vector<IssmDouble>* pf = NULL;
20  Vector<IssmDouble>* df = NULL;
21  Vector<IssmDouble>* ys = NULL;
22 
23  /*parameters:*/
24  int max_nonlinear_iterations;
25  IssmDouble eps_res,eps_rel,eps_abs;
26  HydrologyGlaDSAnalysis* analysis = new HydrologyGlaDSAnalysis();
27 
28  /*Recover parameters (FIXME: from Stress balance for now :( )*/
29  femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
30  femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum);
31  femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
32  femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
33  femmodel->UpdateConstraintsx();
34 
35  int count_out=0;
36  bool converged_out=false;
37  int count_in=0;
38  bool converged_in=false;
39 
40  /*Start non-linear iteration using input velocity: */
41  GetSolutionFromInputsx(&ug,femmodel);
42  Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
43 
44  while(!converged_out){
45 
46  count_in=0;
47  converged_in=false;
48 
49  while(!converged_in){
50  /*save pointer to old solution*/
51  delete old_uf;old_uf=uf;
52  delete ug;
53 
54  SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
55  CreateNodalConstraintsx(&ys,femmodel->nodes);
56  Reduceloadx(pf, Kfs, ys); delete Kfs;
57  femmodel->profiler->Start(SOLVER);
58  Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters);
59  femmodel->profiler->Stop(SOLVER);
60  Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
61 
62  convergence(&converged_in,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs); delete Kff; delete pf; delete df;
63  InputUpdateFromSolutionx(femmodel,ug);
64 
65  /*Increase count: */
66  count_in++;
67  if(count_in>=max_nonlinear_iterations && !converged_in){
68  _printf0_(" maximum number of nonlinear iterations of inner loop (" << max_nonlinear_iterations << ") exceeded\n");
69  converged_in = true;
70  }
71  }
72  if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner loop converged in "<<count_in<<" iterations\n");
73 
74  if(VerboseConvergence()) _printf0_(" updating sheet thickness and channels cross section\n");
75  analysis->UpdateSheetThickness(femmodel);
76  analysis->UpdateChannelCrossSection(femmodel);
77 
78  /*Converged if inner loop converged in one solution*/
79  if(count_in==1) converged_out = true;
80 
81  /*Increase count: */
82  count_out++;
83  if(count_out>=max_nonlinear_iterations && !converged_out){
84  _printf0_(" maximum number of nonlinear iterations of outer loop (" << max_nonlinear_iterations << ") exceeded\n");
85  converged_out = true;
86  }
87  }
88 
89  if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count_out<<"x"<<count_in<<"\n");
90 
91  /*clean-up*/
92  delete uf;
93  delete ug;
94  delete old_uf;
95  delete analysis;
96 }

◆ solutionsequence_nonlinear()

void solutionsequence_nonlinear ( FemModel *  femmodel,
bool  conserve_loads 
)

Definition at line 11 of file solutionsequence_nonlinear.cpp.

11  {
12 
13  /*intermediary: */
14  Matrix<IssmDouble>* Kff = NULL;
15  Matrix<IssmDouble>* Kfs = NULL;
16  Vector<IssmDouble>* ug = NULL;
17  Vector<IssmDouble>* uf = NULL;
18  Vector<IssmDouble>* old_uf = NULL;
19  Vector<IssmDouble>* pf = NULL;
20  Vector<IssmDouble>* df = NULL;
21  Vector<IssmDouble>* ys = NULL;
22 
23  Loads* savedloads=NULL;
24  bool converged;
25  int constraints_converged;
26  int num_unstable_constraints;
27  int count;
28 
29  /*parameters:*/
30  int min_mechanical_constraints;
31  int max_nonlinear_iterations;
32  int configuration_type;
33  IssmDouble eps_res,eps_rel,eps_abs;
34 
35  /*Recover parameters: */
36  femmodel->parameters->FindParam(&min_mechanical_constraints,StressbalanceRiftPenaltyThresholdEnum);
37  femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
38  femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum);
39  femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
40  femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
41  femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
42  femmodel->UpdateConstraintsx();
43 
44  /*Were loads requested as output? : */
45  if(conserve_loads){
46  savedloads = static_cast<Loads*>(femmodel->loads->Copy());
47  }
48 
49  count=0;
50  converged=false;
51 
52  /*Start non-linear iteration using input velocity: */
53  GetSolutionFromInputsx(&ug,femmodel);
54  Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
55 
56  /*Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)*/
57  InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
58  InputUpdateFromSolutionx(femmodel,ug);
59 
60  /*allocate the matrices once and reuse them per iteration*/
61  if(femmodel->loads->numrifts == 0){
62  AllocateSystemMatricesx(&Kff,&Kfs,&df,&pf,femmodel);
63  }
64 
65  for(;;){
66 
67  //save pointer to old velocity
68  delete old_uf;old_uf=uf;
69  delete ug;
70 
71  if(femmodel->loads->numrifts){
72  AllocateSystemMatricesx(&Kff,&Kfs,&df,&pf,femmodel);
73  }
74  SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel, true);
75  CreateNodalConstraintsx(&ys,femmodel->nodes);
76  Reduceloadx(pf, Kfs, ys);
77  femmodel->profiler->Start(SOLVER);
78  Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters);
79  femmodel->profiler->Stop(SOLVER);
80 
81  Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
82 
83  convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs);
84  InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
85  InputUpdateFromSolutionx(femmodel,ug);
86 
87  /*Clean up if rifts*/
88  if(femmodel->loads->numrifts){
89  delete Kfs; delete Kff; delete pf; delete df;
90  }
91 
92  ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
93  if(VerboseConvergence()) _printf0_(" number of unstable constraints: " << num_unstable_constraints << "\n");
94 
95  //rift convergence
96  if (!constraints_converged) {
97  if (converged){
98  if (num_unstable_constraints <= min_mechanical_constraints) converged=true;
99  else converged=false;
100  }
101  }
102 
103  /*Increase count: */
104  count++;
105  if(converged==true){
106  femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,StressbalanceConvergenceNumStepsEnum,count));
107  break;
108  }
109  if(count>=max_nonlinear_iterations){
110  _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");
111  converged=true;
112  femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,StressbalanceConvergenceNumStepsEnum,max_nonlinear_iterations));
113  InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
114  InputUpdateFromSolutionx(femmodel,ug);
115  break;
116  }
117 
118  /*Set the matrix entries to zero if we do an other iteration*/
119  if(femmodel->loads->numrifts==0){
120  Kff->SetZero();
121  Kfs->SetZero();
122  df->Set(0);
123  pf->Set(0);
124  }
125  }
126 
127  /*delete matrices after the iteration loop*/
128  if(femmodel->loads->numrifts==0){
129  delete Kff; delete pf; delete df; delete Kfs;
130  }
131 
132  if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count << "\n");
133 
134  /*clean-up*/
135  if(conserve_loads){
136  delete femmodel->loads;
137  int index=femmodel->AnalysisIndex(configuration_type);
138  femmodel->loads_list[index]=savedloads;
139  femmodel->loads=savedloads;
140  }
141  delete uf;
142  delete ug;
143  delete old_uf;
144 }

◆ solutionsequence_newton()

void solutionsequence_newton ( FemModel *  femmodel)

Definition at line 11 of file solutionsequence_newton.cpp.

11  {
12 
13  /*intermediary: */
14  bool converged;
15  int count,newton;
16  IssmDouble kmax;
17  Matrix<IssmDouble>* Kff = NULL;
18  Matrix<IssmDouble>* Kfs = NULL;
19  Matrix<IssmDouble>* Jff = NULL;
20  Vector<IssmDouble>* ug = NULL;
21  Vector<IssmDouble>* old_ug = NULL;
22  Vector<IssmDouble>* uf = NULL;
23  Vector<IssmDouble>* old_uf = NULL;
24  Vector<IssmDouble>* duf = NULL;
25  Vector<IssmDouble>* pf = NULL;
26  Vector<IssmDouble>* pJf = NULL;
27  Vector<IssmDouble>* df = NULL;
28  Vector<IssmDouble>* ys = NULL;
29 
30  /*parameters:*/
31  int max_nonlinear_iterations;
32  IssmDouble eps_res,eps_rel,eps_abs;
33 
34  /*Recover parameters: */
35  femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
36  femmodel->parameters->FindParam(&newton,StressbalanceIsnewtonEnum);
37  femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum);
38  femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
39  femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
40  femmodel->UpdateConstraintsx();
41 
42  count=1;
43  converged=false;
44 
45  /*Start non-linear iteration using input velocity: */
46  GetSolutionFromInputsx(&ug,femmodel);
47  Reducevectorgtofx(&uf,ug,femmodel->nodes,femmodel->parameters);
48 
49  //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
50  InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
51  InputUpdateFromSolutionx(femmodel,ug);
52 
53  for(;;){
54 
55  delete old_ug;old_ug=ug;
56  delete old_uf;old_uf=uf;
57 
58  /*Solver forward model*/
59  if(count==1 || newton==2){
60  SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
61  CreateNodalConstraintsx(&ys,femmodel->nodes);
62  Reduceloadx(pf,Kfs,ys);delete Kfs;
63  femmodel->profiler->Start(SOLVER);
64  Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);delete df; delete Kff; delete pf;
65  femmodel->profiler->Stop(SOLVER);
66  Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
67  InputUpdateFromSolutionx(femmodel,ug);
68  delete old_ug;old_ug=ug;
69  delete old_uf;old_uf=uf;
70  }
71  uf=old_uf->Duplicate(); old_uf->Copy(uf);
72 
73  /*Prepare next iteration using Newton's method*/
74  SystemMatricesx(&Kff,&Kfs,&pf,&df,&kmax,femmodel);delete df;
75  CreateNodalConstraintsx(&ys,femmodel->nodes);
76  Reduceloadx(pf,Kfs,ys);delete Kfs;
77 
78  pJf=pf->Duplicate();
79  Kff->MatMult(uf,pJf);
80  pJf->Scale(-1.0); pJf->AXPY(pf,+1.0);
81 
82  CreateJacobianMatrixx(&Jff,femmodel,kmax);
83  femmodel->profiler->Start(SOLVER);
84  Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); delete Jff; delete pJf;
85  femmodel->profiler->Stop(SOLVER);
86  uf->AXPY(duf, 1.0); delete duf;
87  Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
88  InputUpdateFromSolutionx(femmodel,ug);
89 
90  /*Check convergence*/
91  convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs);
92  delete Kff; delete pf;
93  if(converged==true){
94  break;
95  }
96  if(count>=max_nonlinear_iterations){
97  _printf0_(" maximum number of Newton iterations (" << max_nonlinear_iterations << ") exceeded\n");
98  break;
99  }
100 
101  count++;
102  }
103 
104  if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n");
105 
106  /*clean-up*/
107  delete uf;
108  delete ug;
109  delete old_ug;
110  delete old_uf;
111 }

◆ solutionsequence_fct()

void solutionsequence_fct ( FemModel *  femmodel)

Definition at line 388 of file solutionsequence_fct.cpp.

388  {/*{{{*/
389 
390  /*intermediary: */
391  IssmDouble theta,deltat,dmax;
392  int configuration_type,analysis_type;
393  Matrix<IssmDouble>* K = NULL;
394  Matrix<IssmDouble>* Mc = NULL;
395  Vector<IssmDouble>* p = NULL;
396  Vector<IssmDouble>* ug = NULL;
397  Vector<IssmDouble>* uf = NULL;
398  MasstransportAnalysis* manalysis = NULL;
399  DamageEvolutionAnalysis* danalysis = NULL;
400 
401  /*Recover parameters: */
402  femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
403  femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
404  femmodel->UpdateConstraintsx();
405  theta = 0.5; //0.5==Crank-Nicolson, 1==Backward Euler
406 
407  /*Create analysis*/
408  femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
409  switch(analysis_type){
410  case MasstransportAnalysisEnum:
411  manalysis = new MasstransportAnalysis();
412  manalysis->MassMatrix(&Mc,femmodel);
413  manalysis->FctKMatrix(&K,NULL,femmodel);
414  manalysis->FctPVector(&p,femmodel);
415  break;
416  case DamageEvolutionAnalysisEnum:
417  danalysis = new DamageEvolutionAnalysis();
418  danalysis->MassMatrix(&Mc,femmodel);
419  danalysis->FctKMatrix(&K,NULL,femmodel);
420  break;
421  default: _error_("analysis type " << EnumToStringx(analysis_type) << " not supported for FCT\n");
422  }
423  delete manalysis;
424  delete danalysis;
425 
426  #ifdef _HAVE_PETSC_
427 
428  /*Convert matrices to PETSc matrices*/
429  Mat D_petsc = NULL;
430  Mat LHS = NULL;
431  Vec RHS = NULL;
432  Vec u = NULL;
433  Vec udot = NULL;
434  Mat K_petsc = K->pmatrix->matrix;
435  Mat Mc_petsc = Mc->pmatrix->matrix;
436  Vec p_petsc = p->pvector->vector;
437 
438  /*Get previous solution u^n*/
439  GetSolutionFromInputsx(&ug,femmodel);
440  Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
441  delete ug;
442 
443  /*Compute lumped mass matrix*/
444  Vec Ml_petsc = NULL;
445  VecDuplicate(uf->pvector->vector,&Ml_petsc);
446  MatGetRowSum(Mc_petsc,Ml_petsc);
447 
448  /*Create D Matrix*/
449  CreateDMatrix(&D_petsc,K_petsc);
450 
451  /*Create LHS: [ML - theta*detlat *(K+D)^n+1]*/
452  CreateLHS(&LHS,&dmax,K_petsc,D_petsc,Ml_petsc,theta,deltat,femmodel,configuration_type);
453 
454  /*Create RHS: [ML + (1 - theta) deltaT L^n] u^n */
455  CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,p_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
456  delete uf;
457  delete p;
458 
459  /*Go solve lower order solution*/
460  femmodel->profiler->Start(SOLVER);
461  SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters);
462  femmodel->profiler->Stop(SOLVER);
463  MatFree(&LHS);
464  VecFree(&RHS);
465 
466  /*Richardson to calculate udot*/
467  RichardsonUdot(&udot,u,Ml_petsc,K_petsc,Mc_petsc);
468  delete K;
469 
470  /*Serialize u and udot*/
471  IssmDouble* udot_serial = NULL;
472  IssmDouble* u_serial = NULL;
473  IssmDouble* ml_serial = NULL;
474  VecToMPISerial(&udot_serial,udot ,IssmComm::GetComm());
475  VecToMPISerial(&u_serial ,u ,IssmComm::GetComm());
476  VecToMPISerial(&ml_serial ,Ml_petsc,IssmComm::GetComm());
477 
478  /*Anti diffusive fluxes*/
479  Vec Fbar = NULL;
480  IssmDouble *Ri_minus_serial = NULL;
481  IssmDouble *Ri_plus_serial = NULL;
482  IssmDouble *ulmin = NULL;
483  IssmDouble *ulmax = NULL;
484  femmodel->GetInputLocalMinMaxOnNodesx(&ulmin,&ulmax,u_serial);
485  CreateRis(&Ri_plus_serial,&Ri_minus_serial,Mc_petsc,D_petsc,ml_serial,u,u_serial,udot_serial,ulmin,ulmax,deltat);
486  CreateFbar(&Fbar,Ri_plus_serial,Ri_minus_serial,Mc_petsc,D_petsc,udot_serial,u_serial,u);
487  xDelete<IssmDouble>(Ri_plus_serial);
488  xDelete<IssmDouble>(Ri_minus_serial);
489  xDelete<IssmDouble>(ulmin);
490  xDelete<IssmDouble>(ulmax);
491 
492  /*Clean up*/
493  MatFree(&D_petsc);
494  delete Mc;
495  xDelete<IssmDouble>(udot_serial);
496  xDelete<IssmDouble>(u_serial);
497  xDelete<IssmDouble>(ml_serial);
498 
499  /*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
500  UpdateSolution(u,udot,Ml_petsc,Fbar,deltat);
501  uf =new Vector<IssmDouble>(u);
502  VecFree(&u);
503  VecFree(&Fbar);
504  VecFree(&udot);
505  VecFree(&Ml_petsc);
506 
507  /*Update Element inputs*/
508  InputUpdateFromSolutionx(femmodel,uf);
509  delete uf;
510 
511  #else
512  _error_("PETSc needs to be installed");
513  #endif
514 }/*}}}*/

◆ solutionsequence_FScoupling_nonlinear()

void solutionsequence_FScoupling_nonlinear ( FemModel *  femmodel,
bool  conserve_loads 
)

Definition at line 11 of file solutionsequence_stokescoupling_nonlinear.cpp.

11  {
12 
13  /*intermediary: */
14  Matrix<IssmDouble> *Kff_horiz = NULL;
15  Matrix<IssmDouble> *Kfs_horiz = NULL;
16  Vector<IssmDouble> *ug_horiz = NULL;
17  Vector<IssmDouble> *uf_horiz = NULL;
18  Vector<IssmDouble> *old_uf_horiz = NULL;
19  Vector<IssmDouble> *pf_horiz = NULL;
20  Vector<IssmDouble> *df_horiz = NULL;
21  Matrix<IssmDouble> *Kff_vert = NULL;
22  Matrix<IssmDouble> *Kfs_vert = NULL;
23  Vector<IssmDouble> *ug_vert = NULL;
24  Vector<IssmDouble> *uf_vert = NULL;
25  Vector<IssmDouble> *pf_vert = NULL;
26  Vector<IssmDouble> *df_vert = NULL;
27  Vector<IssmDouble> *ys = NULL;
28  bool converged;
29  int count;
30 
31  /*parameters:*/
32  int min_mechanical_constraints;
33  int max_nonlinear_iterations;
34  IssmDouble eps_res,eps_rel,eps_abs;
35 
36  /*Recover parameters: */
37  femmodel->parameters->FindParam(&min_mechanical_constraints,StressbalanceRiftPenaltyThresholdEnum);
38  femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
39  femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum);
40  femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
41  femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
42  femmodel->UpdateConstraintsx();
43 
44  count=1;
45  converged=false;
46 
47  /*First get ug_horiz:*/
48  femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
49  GetSolutionFromInputsx(&ug_horiz,femmodel);
50  Reducevectorgtofx(&uf_horiz, ug_horiz, femmodel->nodes,femmodel->parameters);
51 
52  for(;;){
53 
54  /*First stressbalance horiz:*/
55  femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
56 
57  //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
58  InputUpdateFromSolutionx(femmodel,ug_horiz);
59  delete ug_horiz;
60 
61  //save pointer to old velocity
62  delete old_uf_horiz; old_uf_horiz=uf_horiz;
63 
64  /*solve: */
65  SystemMatricesx(&Kff_horiz, &Kfs_horiz, &pf_horiz, &df_horiz, NULL,femmodel);
66  CreateNodalConstraintsx(&ys,femmodel->nodes);
67  Reduceloadx(pf_horiz, Kfs_horiz, ys); delete Kfs_horiz;
68  femmodel->profiler->Start(SOLVER);
69  Solverx(&uf_horiz, Kff_horiz, pf_horiz, old_uf_horiz, df_horiz,femmodel->parameters);
70  femmodel->profiler->Stop(SOLVER);
71  Mergesolutionfromftogx(&ug_horiz, uf_horiz,ys,femmodel->nodes,femmodel->parameters); delete ys;
72  InputUpdateFromSolutionx(femmodel,ug_horiz);
73 
74  convergence(&converged,Kff_horiz,pf_horiz,uf_horiz,old_uf_horiz,eps_res,eps_rel,eps_abs); delete Kff_horiz; delete pf_horiz; delete df_horiz;
75 
76  /*Second compute vertical velocity: */
77  femmodel->SetCurrentConfiguration(StressbalanceVerticalAnalysisEnum);
78 
79  /*solve: */
80  SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert, &df_vert,NULL,femmodel);
81  CreateNodalConstraintsx(&ys,femmodel->nodes);
82  Reduceloadx(pf_vert, Kfs_vert, ys); delete Kfs_vert;
83  femmodel->profiler->Start(SOLVER);
84  Solverx(&uf_vert, Kff_vert, pf_vert, NULL, df_vert,femmodel->parameters); delete Kff_vert; delete pf_vert; delete df_vert;
85  femmodel->profiler->Stop(SOLVER);
86  Mergesolutionfromftogx(&ug_vert, uf_vert,ys,femmodel->nodes,femmodel->parameters);
87  delete uf_vert;
88  delete ys;
89  InputUpdateFromSolutionx(femmodel,ug_vert);
90  delete ug_vert;
91 
92  /*Increase count: */
93  count++;
94  if(converged==true)break;
95  if(count>=max_nonlinear_iterations){
96  _printf0_(" maximum number of iterations (" << max_nonlinear_iterations << ") exceeded\n");
97  break;
98  }
99  }
100 
101  /*clean-up*/
102  delete old_uf_horiz;
103  delete uf_horiz;
104  delete ug_horiz;
105 }

◆ solutionsequence_linear()

void solutionsequence_linear ( FemModel *  femmodel)

Definition at line 10 of file solutionsequence_linear.cpp.

10  {
11 
12  /*intermediary: */
13  Matrix<IssmDouble>* Kff = NULL;
14  Matrix<IssmDouble>* Kfs = NULL;
15  Vector<IssmDouble>* ug = NULL;
16  Vector<IssmDouble>* uf = NULL;
17  Vector<IssmDouble>* pf = NULL;
18  Vector<IssmDouble>* df = NULL;
19  Vector<IssmDouble>* ys = NULL;
20 
21  /*Recover parameters: */
22  femmodel->UpdateConstraintsx();
23  SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
24  CreateNodalConstraintsx(&ys,femmodel->nodes);
25  Reduceloadx(pf, Kfs, ys); delete Kfs;
26 
27  femmodel->profiler->Start(SOLVER);
28  Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters);
29  femmodel->profiler->Stop(SOLVER);
30 
31  /*Clean up*/
32  delete Kff; delete pf; delete df;
33 
34  //#ifdef _HAVE_ADOLC_
35  // for (int i =0; i<uf->svector->M; ++i) {
36  // ADOLC_DUMP_MACRO(uf->svector->vector[i]);
37  // }
38  //#endif
39 
40  Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
41  InputUpdateFromSolutionx(femmodel,ug);
42  delete ug;
43 }

◆ solutionsequence_la()

void solutionsequence_la ( FemModel *  femmodel)

Definition at line 10 of file solutionsequence_la.cpp.

10  {
11 
12  /*intermediary: */
13  Matrix<IssmDouble>* Kff = NULL;
14  Matrix<IssmDouble>* Kfs = NULL;
15  Vector<IssmDouble>* ug = NULL;
16  Vector<IssmDouble>* uf = NULL;
17  Vector<IssmDouble>* pf = NULL;
18  Vector<IssmDouble>* df = NULL;
19  Vector<IssmDouble>* ys = NULL;
20  Vector<IssmDouble>* pug = NULL;
21  Vector<IssmDouble>* pug_old = NULL;
22  IssmDouble eps_rel,r,theta; // 0<theta<.5 -> .15<theta<.45
23  int max_nonlinear_iterations;
24  bool vel_converged = false;
25  bool pressure_converged = false;
26 
27  /*Create analysis*/
28  StressbalanceAnalysis* stressanalysis = new StressbalanceAnalysis();
29  UzawaPressureAnalysis* pressureanalysis = new UzawaPressureAnalysis();
30  femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
31 
32  /*Recover parameters: */
33  femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
34  femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
35 
36  /*Update constraints*/
37  femmodel->UpdateConstraintsx();
38 
39  /*Convergence criterion*/
40  int count = 0;
41  int count_local = 0;
42  Vector<IssmDouble>* vel = NULL;
43  Vector<IssmDouble>* vel_old = NULL;
44  Vector<IssmDouble>* vel_old_local = NULL;
45  GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexPIdEnum);
46  GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexPIdEnum);
47 
48  while(true){
49  count++;
50 
51  /*save pointer to old velocity*/
52  delete vel_old;vel_old=vel->Duplicate(); vel->Copy(vel_old);
53  delete pug_old;pug_old=pug;
54  pressure_converged=false; vel_converged=false;
55 
56  while(true){
57  count_local++;
58  delete vel_old_local;vel_old_local=vel;
59  /*Solve KU=F*/
60  femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
61  SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
62  CreateNodalConstraintsx(&ys,femmodel->nodes);
63  Reduceloadx(pf, Kfs, ys); delete Kfs;
64 
65  femmodel->profiler->Start(SOLVER);
66  Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters);
67  femmodel->profiler->Stop(SOLVER);
68 
69  delete Kff; delete pf; delete df;
70  Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
71 
72  /*Update solution*/
73  InputUpdateFromSolutionx(femmodel,ug); delete ug;
74  GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexPIdEnum);
75  /*Check for convergence*/
76  Vector<IssmDouble>* dvel_local=vel_old_local->Duplicate(); vel_old_local->Copy(dvel_local); dvel_local->AYPX(vel,-1.0);
77  IssmDouble ndu=dvel_local->Norm(NORM_TWO); delete dvel_local;
78  IssmDouble nu =vel_old_local->Norm(NORM_TWO);
79  if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
80  if((ndu/nu)<eps_rel){
81  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " < " << eps_rel*100 << " %\n");
82  break;
83  }
84  else{
85  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << eps_rel*100 << " %\n");
86  }
87  if(count_local>=max_nonlinear_iterations){
88  _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");
89  break;
90  }
91  }
92  count_local=0;
93 
94  femmodel->SetCurrentConfiguration(UzawaPressureAnalysisEnum);
95  SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
96  CreateNodalConstraintsx(&ys,femmodel->nodes);
97  Reduceloadx(pf, Kfs, ys); delete Kfs;
98 
99  femmodel->profiler->Start(SOLVER);
100  Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters);
101  femmodel->profiler->Stop(SOLVER);
102 
103  delete Kff; delete pf; delete df;
104  Mergesolutionfromftogx(&pug, uf,ys,femmodel->nodes,femmodel->parameters); delete uf; delete ys;
105 
106  /*Update solution*/
107  InputUpdateFromSolutionx(femmodel,pug); delete pug;
108  GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexPIdEnum);
109 
110  /*Check for convergence*/
111  Vector<IssmDouble>* dvel=vel_old_local->Duplicate(); vel_old_local->Copy(dvel); dvel->AYPX(vel,-1.0);
112  IssmDouble ndu=dvel->Norm(NORM_TWO); delete dvel;
113  IssmDouble nu =vel_old_local->Norm(NORM_TWO);
114  if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
115  if((ndu/nu)<eps_rel){
116  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " < " << eps_rel*100 << " %\n");
117  vel_converged=true;
118  }
119  else{
120  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << eps_rel*100 << " %\n");
121  vel_converged=false;
122  }
123  Vector<IssmDouble>* dup=pug_old->Duplicate(); pug_old->Copy(dup); dup->AYPX(pug,-1.0);
124  IssmDouble ndp=dup->Norm(NORM_TWO); delete dup;
125  IssmDouble np =pug_old->Norm(NORM_TWO);
126  if (xIsNan<IssmDouble>(ndp) || xIsNan<IssmDouble>(np)) _error_("convergence criterion is NaN!");
127  if((ndp/np)<eps_rel){
128  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: norm(dp)/norm(p)" << ndp/np*100 << " < " << eps_rel*100 << " %\n");
129  pressure_converged=true;
130  }
131  else{
132  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: norm(dp/)/norm(p)" << ndp/np*100 << " > " << eps_rel*100 << " %\n");
133  pressure_converged=false;
134  }
135 
136  if(pressure_converged && vel_converged) break;
137  if(count>=max_nonlinear_iterations){
138  _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");
139  break;
140  }
141  }
142 
143  if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n");
144 
145  delete pug;
146  delete pug_old;
147  delete vel;
148  delete vel_old;
149  delete vel_old_local;
150  delete stressanalysis;
151  delete pressureanalysis;
152 }

◆ solutionsequence_la_theta()

void solutionsequence_la_theta ( FemModel *  femmodel)

Definition at line 10 of file solutionsequence_la_theta.cpp.

10  {
11 
12  /*intermediary: */
13  Matrix<IssmDouble>* Kff = NULL;
14  Matrix<IssmDouble>* Kfs = NULL;
15  Vector<IssmDouble>* ug_old = NULL;
16  Vector<IssmDouble>* ug = NULL;
17  Vector<IssmDouble>* uf = NULL;
18  Vector<IssmDouble>* pf = NULL;
19  Vector<IssmDouble>* df = NULL;
20  Vector<IssmDouble>* ys = NULL;
21  IssmDouble eps_rel,r,theta; // 0<theta<.5 -> .15<theta<.45
22  int max_nonlinear_iterations;
23 
24  /*Create analysis*/
25  StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
26 
27  /*Recover parameters: */
28  femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
29  femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
30  femmodel->parameters->FindParam(&r,AugmentedLagrangianREnum);
31  femmodel->parameters->FindParam(&theta,AugmentedLagrangianThetaEnum);
32 
33  /*Update constraints and initialize d and tau if necessary*/
34  femmodel->UpdateConstraintsx();
35  analysis->InitializeXTH(femmodel->elements,femmodel->parameters);
36 
37  /*Convergence criterion*/
38  int count = 0;
39  GetSolutionFromInputsx(&ug,femmodel);
40  Vector<IssmDouble>* vx = NULL;
41  Vector<IssmDouble>* vx_old = NULL;
42  GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexPIdEnum);
43 
44  while(true){
45  count++;
46 
47  /*save pointer to old velocity*/
48  delete ug_old;ug_old=ug;
49  delete vx_old;vx_old=vx;
50 
51  /*Calculate d*/
52  if(theta>0.){
53  analysis->InputUpdateFromSolutionFSXTH_d(femmodel->elements,femmodel->parameters);
54  }
55 
56  /*Solve KU=F*/
57  SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
58  CreateNodalConstraintsx(&ys,femmodel->nodes);
59  Reduceloadx(pf, Kfs, ys); delete Kfs;
60 
61  femmodel->profiler->Start(SOLVER);
62  Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters);
63  femmodel->profiler->Stop(SOLVER);
64 
65  delete Kff; delete pf; delete df;
66  Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
67 
68  /*Update solution*/
69  InputUpdateFromSolutionx(femmodel,ug);
70 
71  /*Update d and tau accordingly*/
72  analysis->InputUpdateFromSolutionFSXTH_d( femmodel->elements,femmodel->parameters);
73  analysis->InputUpdateFromSolutionFSXTH_tau(femmodel->elements,femmodel->parameters);
74  GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexPIdEnum);
75 
76  /*Check for convergence*/
77  //Vector<IssmDouble>* dug=ug_old->Duplicate(); ug_old->Copy(dug); dug->AYPX(ug,-1.0);
78  //IssmDouble ndu=dug->Norm(NORM_TWO); delete dug;
79  //IssmDouble nu =ug_old->Norm(NORM_TWO);
80  Vector<IssmDouble>* dvx=vx_old->Duplicate(); vx_old->Copy(dvx); dvx->AYPX(vx,-1.0);
81  IssmDouble ndu=dvx->Norm(NORM_TWO); delete dvx;
82  IssmDouble nu =vx_old->Norm(NORM_TWO);
83  if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
84  if((ndu/nu)<eps_rel){
85  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " < " << eps_rel*100 << " %\n");
86  break;
87  }
88  else{
89  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << eps_rel*100 << " %\n");
90  }
91 
92  if(count>=max_nonlinear_iterations){
93  _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");
94  break;
95  }
96  }
97 
98  if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n");
99 
100  delete ug;
101  delete ug_old;
102  delete vx;
103  delete vx_old;
104  delete analysis;
105 }

◆ solutionsequence_adjoint_linear()

void solutionsequence_adjoint_linear ( FemModel *  femmodel)

Definition at line 10 of file solutionsequence_adjoint_linear.cpp.

10  {
11  /*This is axactly the same solutionsequence as solutionsequence_linear except that Reduceloadfromgtofx and Mergesolutionfromftogx
12  * use the flag "true" so that all spc are taken as 0*/
13 
14  /*intermediary: */
15  Matrix<IssmDouble>* Kff = NULL;
16  Matrix<IssmDouble>* Kfs = NULL;
17  Vector<IssmDouble>* ug = NULL;
18  Vector<IssmDouble>* uf = NULL;
19  Vector<IssmDouble>* pf = NULL;
20  Vector<IssmDouble>* df = NULL;
21  Vector<IssmDouble>* ys = NULL;
22 
23  /*Recover parameters: */
24  femmodel->UpdateConstraintsx();
25 
26  SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL,femmodel);
27  CreateNodalConstraintsx(&ys,femmodel->nodes);
28  Reduceloadx(pf, Kfs, ys,true); delete Kfs; //true means spc = 0
29 
30  femmodel->profiler->Start(SOLVER);
31  Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); delete Kff; delete pf; delete df;
32  femmodel->profiler->Stop(SOLVER);
33 
34  Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters,true); delete ys; //true means spc0
35  InputUpdateFromSolutionx(femmodel,ug);
36  delete ug; delete uf;
37 }

◆ solutionsequence_schurcg()

void solutionsequence_schurcg ( FemModel *  femmodel)

Definition at line 834 of file solutionsequence_schurcg.cpp.

834 {_error_("PETSc needs to be installed");}

◆ convergence()

void convergence ( bool *  pconverged,
Matrix< IssmDouble > *  K_ff,
Vector< IssmDouble > *  p_f,
Vector< IssmDouble > *  u_f,
Vector< IssmDouble > *  u_f_old,
IssmDouble  eps_res,
IssmDouble  eps_rel,
IssmDouble  eps_abs 
)

Definition at line 9 of file convergence.cpp.

9  {
10 
11  /*output*/
12  bool converged=false;
13 
14  /*intermediary*/
15  Vector<IssmDouble>* KUold=NULL;
16  Vector<IssmDouble>* KUoldF=NULL;
17  Vector<IssmDouble>* duf=NULL;
18  IssmDouble ndu,nduinf,nu;
19  IssmDouble nKUoldF;
20  IssmDouble nF;
21  IssmDouble res;
22  int analysis_type;
23 
24  if(VerboseModule()) _printf0_(" checking convergence\n");
25 
26  /*If uf is NULL in input, f-set is nil, model is fully constrained, therefore converged from
27  * the get go: */
28  if(uf->IsEmpty()){
29  *pconverged=true;
30  return;
31  }
32 
33  /*Force equilibrium (Mandatory)*/
34 
35  /*compute K[n]U[n-1]F = K[n]U[n-1] - F*/
36  _assert_(uf); _assert_(Kff);
37  KUold=uf->Duplicate(); Kff->MatMult(old_uf,KUold);
38  KUoldF=KUold->Duplicate();KUold->Copy(KUoldF); KUoldF->AYPX(pf,-1.0);
39  nKUoldF=KUoldF->Norm(NORM_TWO);
40  nF=pf->Norm(NORM_TWO);
41  res=nKUoldF/nF;
42  if (xIsNan<IssmDouble>(res)){
43  _printf0_("norm nf = " << nF << "f and norm kuold = " << nKUoldF << "f\n");
44  _error_("mechanical equilibrium convergence criterion is NaN!");
45  }
46 
47  //clean up
48  delete KUold;
49  delete KUoldF;
50 
51  //print
52  if(res<eps_res){
53  if(VerboseConvergence()) _printf0_(setw(50)<<left<<" mechanical equilibrium convergence criterion"<<res*100<< " < "<<eps_res*100<<" %\n");
54  converged=true;
55  }
56  else{
57  if(VerboseConvergence()) _printf0_(setw(50)<<left<<" mechanical equilibrium convergence criterion"<<res*100<<" > "<<eps_res*100<<" %\n");
58  converged=false;
59  }
60 
61  /*Relative criterion (optional)*/
62  if (!xIsNan<IssmDouble>(eps_rel) || (VerboseConvergence())){
63 
64  //compute norm(du)/norm(u)
65  duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0);
66  ndu=duf->Norm(NORM_TWO); nu=old_uf->Norm(NORM_TWO);
67 
68  if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
69 
70  //clean up
71  delete duf;
72 
73  //print
74  if (!xIsNan<IssmDouble>(eps_rel)){
75  if((ndu/nu)<eps_rel){
76  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " < " << eps_rel*100 << " %\n");
77  }
78  else{
79  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << eps_rel*100 << " %\n");
80  converged=false;
81  }
82  }
83  else _printf0_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " %\n");
84 
85  }
86 
87  /*Absolute criterion (Optional) = max(du)*/
88  if (!xIsNan<IssmDouble>(eps_abs) || (VerboseConvergence())){
89 
90  //compute max(du)
91  duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0);
92  ndu=duf->Norm(NORM_TWO); nduinf=duf->Norm(NORM_INF);
93  if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
94 
95  //clean up
96  delete duf;
97 
98  //print
99  if (!xIsNan<IssmDouble>(eps_abs)){
100  if ((nduinf)<eps_abs){
101  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: max(du)" << nduinf << " < " << eps_abs << "\n");
102  }
103  else{
104  if(VerboseConvergence()) _printf0_(setw(50) << left << " Convergence criterion: max(du)" << nduinf << " > " << eps_abs << "\n");
105  converged=false;
106  }
107  }
108  else _printf0_(setw(50) << left << " Convergence criterion: max(du)" << nduinf << "\n");
109 
110  }
111 
112  /*assign output*/
113  *pconverged=converged;
114 }
DataSet::Size
int Size()
Definition: DataSet.cpp:399
Matrix< IssmDouble >
SOLVER
#define SOLVER
Definition: Profiler.h:16
FemModel::AnalysisIndex
int AnalysisIndex(int)
Definition: FemModel.cpp:220
StressbalanceAnalysis
Definition: StressbalanceAnalysis.h:11
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
EnthalpyAnalysis::UpdateBasalConstraints
static void UpdateBasalConstraints(FemModel *femmodel)
Definition: EnthalpyAnalysis.cpp:1670
MasstransportAnalysisEnum
@ MasstransportAnalysisEnum
Definition: EnumDefinitions.h:1163
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
StressbalanceAnalysisEnum
@ StressbalanceAnalysisEnum
Definition: EnumDefinitions.h:1285
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
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
DamageEvolutionAnalysisEnum
@ DamageEvolutionAnalysisEnum
Definition: EnumDefinitions.h:1026
ConvergedEnum
@ ConvergedEnum
Definition: EnumDefinitions.h:514
ThermalIsdynamicbasalspcEnum
@ ThermalIsdynamicbasalspcEnum
Definition: EnumDefinitions.h:416
HydrologydcIsefficientlayerEnum
@ HydrologydcIsefficientlayerEnum
Definition: EnumDefinitions.h:185
GetVectorFromInputsx
void GetVectorFromInputsx(IssmDouble **pvector, int *pvector_size, FemModel *femmodel, int name)
Definition: GetVectorFromInputsx.cpp:81
FemModel::results
Results * results
Definition: FemModel.h:48
StressbalanceConvergenceNumStepsEnum
@ StressbalanceConvergenceNumStepsEnum
Definition: EnumDefinitions.h:1286
HydrologyGlaDSAnalysis::UpdateSheetThickness
void UpdateSheetThickness(FemModel *femmodel)
Definition: HydrologyGlaDSAnalysis.cpp:458
NORM_INF
@ NORM_INF
Definition: toolkitsenums.h:15
EnthalpyAnalysis::ComputeBasalMeltingrate
static void ComputeBasalMeltingrate(FemModel *femmodel)
Definition: EnthalpyAnalysis.cpp:355
PressureEnum
@ PressureEnum
Definition: EnumDefinitions.h:664
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
Matrix::MatMult
void MatMult(Vector< doubletype > *X, Vector< doubletype > *AX)
Definition: Matrix.h:239
AllocateSystemMatricesx
void AllocateSystemMatricesx(Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, Vector< IssmDouble > **pdf, Vector< IssmDouble > **ppf, FemModel *femmodel)
Definition: AllocateSystemMatricesx.cpp:9
ThermalPenaltyThresholdEnum
@ ThermalPenaltyThresholdEnum
Definition: EnumDefinitions.h:422
HydrologydcRelTolEnum
@ HydrologydcRelTolEnum
Definition: EnumDefinitions.h:190
StressbalanceRiftPenaltyThresholdEnum
@ StressbalanceRiftPenaltyThresholdEnum
Definition: EnumDefinitions.h:413
VerboseModule
bool VerboseModule(void)
Definition: Verbosity.cpp:23
AugmentedLagrangianThetaEnum
@ AugmentedLagrangianThetaEnum
Definition: EnumDefinitions.h:41
solutionsequence_linear
void solutionsequence_linear(FemModel *femmodel)
Definition: solutionsequence_linear.cpp:10
NORM_TWO
@ NORM_TWO
Definition: toolkitsenums.h:15
VecFree
void VecFree(Vec *pvec)
Definition: VecFree.cpp:16
HydrologyDCEfficientAnalysis
Definition: HydrologyDCEfficientAnalysis.h:12
StressbalanceMaxiterEnum
@ StressbalanceMaxiterEnum
Definition: EnumDefinitions.h:407
SolverxPetsc
void SolverxPetsc(Vec *puf, Mat Kff, Vec pf, Vec uf0, Vec df, Parameters *parameters)
Definition: PetscSolver.cpp:35
Loads::Copy
Loads * Copy()
Definition: Loads.cpp:42
HydrologyDCInefficientAnalysisEnum
@ HydrologyDCInefficientAnalysisEnum
Definition: EnumDefinitions.h:1099
Vector::Norm
doubletype Norm(NormMode norm_type)
Definition: Vector.h:342
Matrix::SetZero
void SetZero(void)
Definition: Matrix.h:312
Vector::Scale
void Scale(doubletype scale_factor)
Definition: Vector.h:355
HydrologyDCEfficientAnalysis::InitZigZagCounter
void InitZigZagCounter(FemModel *femmodel)
Definition: HydrologyDCEfficientAnalysis.cpp:136
UzawaPressureAnalysisEnum
@ UzawaPressureAnalysisEnum
Definition: EnumDefinitions.h:1320
StressbalanceIsnewtonEnum
@ StressbalanceIsnewtonEnum
Definition: EnumDefinitions.h:406
Parameters::SetParam
void SetParam(bool boolean, int enum_type)
Definition: Parameters.cpp:441
Vector::AXPY
void AXPY(Vector *X, doubletype a)
Definition: Vector.h:256
MasstransportAnalysis::MassMatrix
void MassMatrix(Matrix< IssmDouble > **pMff, FemModel *femmodel)
Definition: MasstransportAnalysis.cpp:1176
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
ThermalReltolEnum
@ ThermalReltolEnum
Definition: EnumDefinitions.h:423
HydrologyDCInefficientAnalysis::ElementizeIdsMask
void ElementizeIdsMask(FemModel *femmodel)
Definition: HydrologyDCInefficientAnalysis.cpp:816
VecToMPISerial
int VecToMPISerial(double **pgathered_vector, Vec vector, ISSM_MPI_Comm comm, bool broadcast=true)
Definition: VecToMPISerial.cpp:14
FemModel::GetInputLocalMinMaxOnNodesx
void GetInputLocalMinMaxOnNodesx(IssmDouble **pmin, IssmDouble **pmax, IssmDouble *ug)
Definition: FemModel.cpp:1308
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
StressbalanceAbstolEnum
@ StressbalanceAbstolEnum
Definition: EnumDefinitions.h:404
CreateJacobianMatrixx
void CreateJacobianMatrixx(Matrix< IssmDouble > **pJff, FemModel *femmodel, IssmDouble kmax)
Definition: CreateJacobianMatrixx.cpp:10
ThermalIsenthalpyEnum
@ ThermalIsenthalpyEnum
Definition: EnumDefinitions.h:417
HydrologyDCInefficientAnalysis
Definition: HydrologyDCInefficientAnalysis.h:12
StressbalanceReltolEnum
@ StressbalanceReltolEnum
Definition: EnumDefinitions.h:410
HydrologydcMaxIterEnum
@ HydrologydcMaxIterEnum
Definition: EnumDefinitions.h:187
HydrologyDCEfficientAnalysis::ComputeEPLThickness
void ComputeEPLThickness(FemModel *femmodel)
Definition: HydrologyDCEfficientAnalysis.cpp:615
HydrologyGlaDSAnalysis::UpdateChannelCrossSection
void UpdateChannelCrossSection(FemModel *femmodel)
Definition: HydrologyGlaDSAnalysis.cpp:607
ConfigurationTypeEnum
@ ConfigurationTypeEnum
Definition: EnumDefinitions.h:101
Results::AddResult
int AddResult(ExternalResult *result)
Definition: Results.cpp:33
DamageEvolutionAnalysis::FctKMatrix
void FctKMatrix(Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, FemModel *femmodel)
Definition: DamageEvolutionAnalysis.cpp:868
InputUpdateFromConstantx
void InputUpdateFromConstantx(FemModel *femmodel, bool constant, int name)
Definition: InputUpdateFromConstantx.cpp:10
FemModel::loads
Loads * loads
Definition: FemModel.h:54
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
FemModel::elements
Elements * elements
Definition: FemModel.h:44
Reducevectorgtofx
void Reducevectorgtofx(Vector< IssmDouble > **puf, Vector< IssmDouble > *ug, Nodes *nodes, Parameters *parameters)
Definition: Reducevectorgtofx.cpp:8
StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_tau
void InputUpdateFromSolutionFSXTH_tau(Elements *elements, Parameters *parameters)
Definition: StressbalanceAnalysis.cpp:5535
MasstransportAnalysis::FctKMatrix
void FctKMatrix(Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, FemModel *femmodel)
Definition: MasstransportAnalysis.cpp:1128
VertexPIdEnum
@ VertexPIdEnum
Definition: EnumDefinitions.h:1324
Reduceloadx
void Reduceloadx(Vector< IssmDouble > *pf, Matrix< IssmDouble > *Kfs, Vector< IssmDouble > *y_s, bool flag_ys0)
Definition: Reduceloadx.cpp:14
Loads::numrifts
int numrifts
Definition: Loads.h:20
GenericExternalResult
Definition: GenericExternalResult.h:21
Loads
Declaration of Loads class.
Definition: Loads.h:16
StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_d
void InputUpdateFromSolutionFSXTH_d(Elements *elements, Parameters *parameters)
Definition: StressbalanceAnalysis.cpp:5322
DamageEvolutionAnalysis
Definition: DamageEvolutionAnalysis.h:11
L2ProjectionEPLAnalysisEnum
@ L2ProjectionEPLAnalysisEnum
Definition: EnumDefinitions.h:1137
UzawaPressureAnalysis
Definition: UzawaPressureAnalysis.h:11
Vector::AYPX
void AYPX(Vector *X, doubletype a)
Definition: Vector.h:267
StressbalanceAnalysis::InitializeXTH
void InitializeXTH(Elements *elements, Parameters *parameters)
Definition: StressbalanceAnalysis.cpp:5155
MatFree
void MatFree(Mat *pmat)
Definition: MatFree.cpp:16
StressbalanceVerticalAnalysisEnum
@ StressbalanceVerticalAnalysisEnum
Definition: EnumDefinitions.h:1289
_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
StressbalanceRestolEnum
@ StressbalanceRestolEnum
Definition: EnumDefinitions.h:412
GetSolutionFromInputsx
void GetSolutionFromInputsx(Vector< IssmDouble > **psolution, FemModel *femmodel)
Definition: GetSolutionFromInputsx.cpp:9
VerboseSolution
bool VerboseSolution(void)
Definition: Verbosity.cpp:24
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
HydrologyDCEfficientAnalysis::ResetCounter
void ResetCounter(FemModel *femmodel)
Definition: HydrologyDCEfficientAnalysis.cpp:143
FemModel::SetCurrentConfiguration
void SetCurrentConfiguration(int configuration_type)
Definition: FemModel.cpp:634
HydrologyGlaDSAnalysis
Definition: HydrologyGlaDSAnalysis.h:11
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
ResetConstraintsx
void ResetConstraintsx(FemModel *femmodel)
Definition: ResetConstraintsx.cpp:16
AugmentedLagrangianREnum
@ AugmentedLagrangianREnum
Definition: EnumDefinitions.h:37
InputToL2ProjectEnum
@ InputToL2ProjectEnum
Definition: EnumDefinitions.h:206
AnalysisTypeEnum
@ AnalysisTypeEnum
Definition: EnumDefinitions.h:36
ThermalMaxiterEnum
@ ThermalMaxiterEnum
Definition: EnumDefinitions.h:418
DamageEvolutionAnalysis::MassMatrix
void MassMatrix(Matrix< IssmDouble > **pMff, FemModel *femmodel)
Definition: DamageEvolutionAnalysis.cpp:896
Vector::Copy
void Copy(Vector *to)
Definition: Vector.h:319
HydrologyDCInefficientAnalysis::ElementizeEplMask
void ElementizeEplMask(FemModel *femmodel)
Definition: HydrologyDCInefficientAnalysis.cpp:756
FemModel::loads_list
Loads ** loads_list
Definition: FemModel.h:55
CreateNodalConstraintsx
void CreateNodalConstraintsx(Vector< IssmDouble > **pys, Nodes *nodes)
Definition: CreateNodalConstraintsx.cpp:10
MasstransportAnalysis::FctPVector
void FctPVector(Vector< IssmDouble > **ppf, FemModel *femmodel)
Definition: MasstransportAnalysis.cpp:1156
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
convergence
void convergence(bool *pconverged, Matrix< IssmDouble > *Kff, Vector< IssmDouble > *pf, Vector< IssmDouble > *uf, Vector< IssmDouble > *old_uf, IssmDouble eps_res, IssmDouble eps_rel, IssmDouble eps_abs)
Definition: convergence.cpp:9
FemModel::profiler
Profiler * profiler
Definition: FemModel.h:42
Vector< IssmDouble >
MeltingOffsetEnum
@ MeltingOffsetEnum
Definition: EnumDefinitions.h:269
MasstransportAnalysis
Definition: MasstransportAnalysis.h:11
Vector::Set
void Set(doubletype value)
Definition: Vector.h:245
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:25 for Ice Sheet System Model by   doxygen 1.8.19