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

: core of the steadystate solution More...

#include "./cores.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"
#include "../solutionsequences/solutionsequences.h"

Go to the source code of this file.

Functions

bool steadystateconvergence (Vector< IssmDouble > *tg, Vector< IssmDouble > *tg_old, Vector< IssmDouble > *ug, Vector< IssmDouble > *ug_old, IssmDouble reltol)
 
void steadystate_core (FemModel *femmodel)
 

Detailed Description

: core of the steadystate solution

Definition in file steadystate_core.cpp.

Function Documentation

◆ steadystateconvergence()

bool steadystateconvergence ( Vector< IssmDouble > *  tg,
Vector< IssmDouble > *  tg_old,
Vector< IssmDouble > *  ug,
Vector< IssmDouble > *  ug_old,
IssmDouble  reltol 
)

Definition at line 91 of file steadystate_core.cpp.

91  {//{{{
92 
93  /*Output*/
94  bool converged = true;
95 
96  /*Intermediary*/
97  Vector<IssmDouble>* dug = NULL;
98  Vector<IssmDouble>* dtg = NULL;
99  IssmDouble ndt,nt;
100  IssmDouble ndu,nu;
101 
102  /*compute norm(du)/norm(u)*/
103  dug=ug_old->Duplicate(); ug_old->Copy(dug); dug->AYPX(ug,-1.0);
104  ndu=dug->Norm(NORM_TWO); nu=ug_old->Norm(NORM_TWO);
105  if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
106  if((ndu/nu)<reltol){
107  if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<" Velocity convergence: norm(du)/norm(u)"<<ndu/nu*100<<" < "<<reltol*100<<" %\n");
108  }
109  else{
110  if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<" Velocity convergence: norm(du)/norm(u)"<<ndu/nu*100<<" > "<<reltol*100<<" %\n");
111  converged=false;
112  }
113 
114  /*compute norm(dt)/norm(t)*/
115  dtg=tg_old->Duplicate(); tg_old->Copy(dtg); dtg->AYPX(tg,-1.0);
116  ndt=dtg->Norm(NORM_TWO); nt=tg_old->Norm(NORM_TWO);
117  if (xIsNan<IssmDouble>(ndt) || xIsNan<IssmDouble>(nt)) _error_("convergence criterion is NaN!");
118  if((ndt/nt)<reltol){
119  if(VerboseConvergence()) _printf0_(setw(50)<<left<<" Temperature convergence: norm(dt)/norm(t)"<<ndt/nt*100<<" < "<<reltol*100<<" %\n");
120  }
121  else{
122  if(VerboseConvergence()) _printf0_(setw(50)<<left<<" Temperature convergence: norm(dt)/norm(t)"<<ndt/nt*100<<" > "<<reltol*100<<" %\n");
123  converged=false;
124  }
125 
126  /*clean up and return*/
127  delete dtg;
128  delete dug;
129  return converged;
130 }//}}}

◆ steadystate_core()

void steadystate_core ( FemModel femmodel)

Definition at line 20 of file steadystate_core.cpp.

20  { //{{{
21 
22  /*intermediary: */
23  int step;
24  Vector<IssmDouble>* ug = NULL;
25  Vector<IssmDouble>* ug_old = NULL;
26  Vector<IssmDouble>* tg = NULL;
27  Vector<IssmDouble>* tg_old = NULL;
28 
29  /*parameters: */
30  bool save_results,isenthalpy;
31  int maxiter;
32  IssmDouble reltol;
33  int numoutputs = 0;
34  char** requested_outputs = NULL;
35 
36  /* recover parameters:*/
43  if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SteadystateRequestedOutputsEnum);
44 
45  /*intialize counters: */
46  step=1;
47 
48  for(;;){
49 
50  /* Compute first velocity, then temperature due to high sensitivity of temperature to velocity. */
51  if(VerboseSolution()) _printf0_("\n======================================================\n");
52  if(VerboseSolution()) _printf0_(" computing velocity and temperature for step: " << step << "\n");
53  if(VerboseSolution()) _printf0_("====================================================\n");
54 
55  if(VerboseSolution()) _printf0_("\n -- computing new velocity -- \n\n");
58 
59  if(VerboseSolution()) _printf0_("\n -- computing new temperature --\n\n");
61  if(!isenthalpy)femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);/*Could be MeltingAnalysis...*/
63 
64  if(step>1){
65  if(VerboseSolution()) _printf0_(" checking steadystate convergence\n");
66  if(steadystateconvergence(tg,tg_old,ug,ug_old,reltol)) break;
67  }
68  if(step>maxiter){
69  if(VerboseSolution()) _printf0_(" maximum number steadystate iterations " << maxiter << " reached\n");
70  break;
71  }
72 
73  /*update results and increase counter*/
74  delete tg_old;tg_old=tg;
75  delete ug_old;ug_old=ug;
76  step++;
77  }
78 
79  if(save_results){
80  if(VerboseSolution()) _printf0_(" saving results\n");
81  femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
82  }
83 
84  /*Free ressources:*/
85  delete tg_old;
86  delete ug_old;
87  delete tg;
88  delete ug;
89  if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
90 }//}}}
SaveResultsEnum
@ SaveResultsEnum
Definition: EnumDefinitions.h:302
IssmDouble
double IssmDouble
Definition: types.h:37
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
VerboseConvergence
bool VerboseConvergence(void)
Definition: Verbosity.cpp:26
thermal_core
void thermal_core(FemModel *femmodel)
Definition: thermal_core.cpp:13
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
steadystateconvergence
bool steadystateconvergence(Vector< IssmDouble > *tg, Vector< IssmDouble > *tg_old, Vector< IssmDouble > *ug, Vector< IssmDouble > *ug_old, IssmDouble reltol)
Definition: steadystate_core.cpp:91
FemModel::results
Results * results
Definition: FemModel.h:48
Vector::Duplicate
Vector< doubletype > * Duplicate(void)
Definition: Vector.h:230
SteadystateNumRequestedOutputsEnum
@ SteadystateNumRequestedOutputsEnum
Definition: EnumDefinitions.h:400
SteadystateRequestedOutputsEnum
@ SteadystateRequestedOutputsEnum
Definition: EnumDefinitions.h:402
NORM_TWO
@ NORM_TWO
Definition: toolkitsenums.h:15
Vector::Norm
doubletype Norm(NormMode norm_type)
Definition: Vector.h:342
Parameters::SetParam
void SetParam(bool boolean, int enum_type)
Definition: Parameters.cpp:441
ThermalIsenthalpyEnum
@ ThermalIsenthalpyEnum
Definition: EnumDefinitions.h:417
ThermalAnalysisEnum
@ ThermalAnalysisEnum
Definition: EnumDefinitions.h:1302
Vector::AYPX
void AYPX(Vector *X, doubletype a)
Definition: Vector.h:267
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
GetSolutionFromInputsx
void GetSolutionFromInputsx(Vector< IssmDouble > **psolution, FemModel *femmodel)
Definition: GetSolutionFromInputsx.cpp:9
VerboseSolution
bool VerboseSolution(void)
Definition: Verbosity.cpp:24
FemModel::RequestedOutputsx
void RequestedOutputsx(Results **presults, char **requested_outputs, int numoutputs, bool save_results=true)
Definition: FemModel.cpp:2267
FemModel::SetCurrentConfiguration
void SetCurrentConfiguration(int configuration_type)
Definition: FemModel.cpp:634
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
SteadystateMaxiterEnum
@ SteadystateMaxiterEnum
Definition: EnumDefinitions.h:399
Vector::Copy
void Copy(Vector *to)
Definition: Vector.h:319
SteadystateReltolEnum
@ SteadystateReltolEnum
Definition: EnumDefinitions.h:401
stressbalance_core
void stressbalance_core(FemModel *femmodel)
Definition: stressbalance_core.cpp:13
Vector< IssmDouble >
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16