Ice Sheet System Model  4.18
Code documentation
Functions
convergence.cpp File Reference

: figure out if convergence has been reached More...

#include "../classes/classes.h"
#include "../modules/modules.h"
#include "../shared/shared.h"

Go to the source code of this file.

Functions

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)
 

Detailed Description

: figure out if convergence has been reached

Definition in file convergence.cpp.

Function Documentation

◆ 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 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 }
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
VerboseConvergence
bool VerboseConvergence(void)
Definition: Verbosity.cpp:26
NORM_INF
@ NORM_INF
Definition: toolkitsenums.h: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
VerboseModule
bool VerboseModule(void)
Definition: Verbosity.cpp:23
NORM_TWO
@ NORM_TWO
Definition: toolkitsenums.h:15
Vector::Norm
doubletype Norm(NormMode norm_type)
Definition: Vector.h:342
Vector::IsEmpty
bool IsEmpty(void)
Definition: Vector.h:196
Vector::AYPX
void AYPX(Vector *X, doubletype a)
Definition: Vector.h:267
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Vector::Copy
void Copy(Vector *to)
Definition: Vector.h:319
Vector< IssmDouble >