Ice Sheet System Model  4.18
Code documentation
sealevelchange_core.cpp
Go to the documentation of this file.
1 
5 #include "./cores.h"
6 #include "../toolkits/toolkits.h"
7 #include "../classes/classes.h"
8 #include "../classes/Inputs2/TriaInput2.h"
9 #include "../classes/Inputs2/TransientInput2.h"
10 #include "../classes/Inputs2/DatasetInput2.h"
11 #include "../shared/shared.h"
12 #include "../modules/modules.h"
13 #include "../solutionsequences/solutionsequences.h"
14 
15 
16 /*main cores:*/
17 void sealevelchange_core(FemModel* femmodel){ /*{{{*/
18 
19  /*Start profiler*/
20  femmodel->profiler->Start(SLRCORE);
21 
22  /*Parameters, variables:*/
23  bool save_results;
24  bool isslr=0;
25  bool isgia=0;
26  int solution_type;
27 
28  /*Retrieve parameters:*/
29  femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
30  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
31  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
32 
33  /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
34  if(solution_type==SealevelriseSolutionEnum){
35  isslr=1;
36  isgia=1;
37  }
38 
39  /*Should we be here?:*/
40  if(!isslr)return;
41 
42  /*Verbose: */
43  if(VerboseSolution()) _printf0_(" computing sea level rise\n");
44 
45  /*Run gia core: */
46  if(isgia){
47  #ifdef _HAVE_GIA_
48  gia_core(femmodel);
49  #else
50  _error_("ISSM was not compiled with gia capabilities. Exiting");
51  #endif
52  }
53 
54  /*set SLR configuration: */
55  femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
56 
57  /*run geometry core: */
58  sealevelrise_core_geometry(femmodel);
59 
60  /*Run geodetic:*/
61  grd_core(femmodel);
62 
63  /*Run steric core for sure:*/
64  dynstr_core(femmodel);
65 
66  /*Save results: */
67  if(save_results){
68  int numoutputs = 0;
69  char **requested_outputs = NULL;
70 
71  if(VerboseSolution()) _printf0_(" saving results\n");
72  femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum);
73  femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
74  if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
75  }
76 
77  /*requested dependents: */
78  if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx();
79 
80  /*End profiler*/
81  femmodel->profiler->Stop(SLRCORE);
82 }
83 /*}}}*/
84 void grd_core(FemModel* femmodel){ /*{{{*/
85 
86  /*Gravity rotation deformation core GRD: */
87 
88  /*variables:*/
89  Vector<IssmDouble> *RSLg = NULL;
90  Vector<IssmDouble> *BPg = NULL;
91  Vector<IssmDouble> *RSLg_rate = NULL;
92  Vector<IssmDouble> *RSLg_eustatic = NULL;
93  Vector<IssmDouble> *U_esa = NULL;
94  Vector<IssmDouble> *U_esa_rate = NULL;
95  Vector<IssmDouble> *N_esa = NULL;
96  Vector<IssmDouble> *N_esa_rate = NULL;
97  Vector<IssmDouble> *U_north_esa = NULL;
98  Vector<IssmDouble> *U_east_esa = NULL;
99  Vector<IssmDouble> *N_gia= NULL;
100  Vector<IssmDouble> *U_gia= NULL;
101  Vector<IssmDouble> *N_gia_rate= NULL;
102  Vector<IssmDouble> *U_gia_rate= NULL;
103  SealevelMasks* masks=NULL;
104 
105  /*parameters:*/
106  bool iscoupler;
107  int solution_type;
108  int modelid,earthid;
109  bool istransientmasstransport;
110  int frequency,count;
111  int horiz;
112  int geodetic=0;
113  IssmDouble dt;
114  IssmDouble oceanarea;
115  int bp_compute_fingerprints=0;
116 
117  /*Should we even be here?:*/
118  femmodel->parameters->FindParam(&geodetic,SolidearthSettingsComputesealevelchangeEnum); if(!geodetic)return;
119 
120  /*Verbose: */
121  if(VerboseSolution()) _printf0_(" computing geodetic sea level rise\n");
122 
123  /*retrieve more parameters:*/
124  femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
125  femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
126  femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum);
127  femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
128  femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
129  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
130 
131  if(iscoupler){
132  femmodel->parameters->FindParam(&modelid,ModelIdEnum);
133  femmodel->parameters->FindParam(&earthid,EarthIdEnum);
134  femmodel->parameters->FindParam(&istransientmasstransport,TransientIsmasstransportEnum);
135  }
136  else{
137  /* we are here, we are not running in a coupler, so we will indeed compute SLR,
138  * so make sure we are identified as being the Earth.:*/
139  modelid=1; earthid=1;
140  /* in addition, if we are running solution_type SealevelriseSolutionEnum, make sure we
141  * run, irresepective of the time settings:*/
142  count=frequency;
143  }
144 
145  /*If we are running in coupled mode, the Earth model needs to run its own mass transport (if
146  * not already done by the mass trasnport module. For ice caps, they rely on the transient mass
147  * transport module exclusively:*/
148  if(iscoupler) if(modelid==earthid) if(!istransientmasstransport) EarthMassTransport(femmodel);
149 
150  /*increment counter, or call solution core if count==frequency:*/
151  if (count<frequency){
152  count++; femmodel->parameters->SetParam(count,SealevelriseRunCountEnum);
153  return;
154  }
155 
156  /*call sea-level rise sub cores:*/
157  if(iscoupler){
158  /*transfer cumulated deltathickness forcing from ice caps to earth model: */
159  TransferForcing(femmodel,SealevelriseCumDeltathicknessEnum);
160 
161  /*we have accumulated thicknesses, dump them in deltathcikness: */
162  if(modelid==earthid)InputDuplicatex(femmodel,SealevelriseCumDeltathicknessEnum,SurfaceloadIceThicknessChangeEnum);
163  }
164 
165  /*run cores:*/
166  if(modelid==earthid){
167 
168  /*call masks core: */
169  masks=sealevelrise_core_masks(femmodel);
170 
171  /*call eustatic core (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */
172  RSLg_eustatic=sealevelrise_core_eustatic(femmodel,masks,&oceanarea);
173 
174  /*call non-eustatic core (ocean loading tems - 2nd and 5th terms on the RHS of Farrel and Clark) */
175  RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,oceanarea);
176 
177  /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */
178  sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,masks);
179 
180  /*compute viscosus (GIA) geodetic signatures:*/
181  sealevelrise_core_viscous(&U_gia,&N_gia,femmodel,RSLg);
182 
183  /*compute sea-level rise (low-order spherical harmonics coefficients) diagnostics:*/
184  sealevelrise_diagnostics(femmodel,RSLg);
185 
186  /*recover N_esa = U_esa + RSLg:*/
187  N_esa=U_esa->Duplicate(); U_esa->Copy(N_esa); N_esa->AXPY(RSLg,1);
188 
189  /*if we had bottom pressure loading, remove dynamic sea level from geoid: */
190  femmodel->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
191  if(bp_compute_fingerprints){
192  GetVectorFromInputsx(&BPg,femmodel,DslSeaWaterPressureChangeAtSeaFloor,VertexSIdEnum);
193  N_esa->AXPY(BPg,-1);
194  }
195 
196  /*transform these values into rates (as we only run this once each frequency turn:*/
197  N_esa_rate=N_esa->Duplicate(); N_esa->Copy(N_esa_rate); N_esa_rate->Scale(1/(dt*frequency));
198  U_esa_rate=U_esa->Duplicate(); U_esa->Copy(U_esa_rate); U_esa_rate->Scale(1/(dt*frequency));
199  N_gia_rate=N_gia->Duplicate(); N_gia->Copy(N_gia_rate); N_gia_rate->Scale(1/(dt*frequency));
200  U_gia_rate=U_gia->Duplicate(); U_gia->Copy(U_gia_rate); U_gia_rate->Scale(1/(dt*frequency));
201  RSLg_rate=RSLg->Duplicate(); RSLg->Copy(RSLg_rate); RSLg_rate->Scale(1/(dt*frequency));
202 
203  /*get some results into elements:{{{*/
204  InputUpdateFromVectorx(femmodel,U_esa_rate,SealevelUEsaRateEnum,VertexSIdEnum);
205  InputUpdateFromVectorx(femmodel,N_esa_rate,SealevelNEsaRateEnum,VertexSIdEnum);
206  InputUpdateFromVectorx(femmodel,U_gia_rate,UGiaRateEnum,VertexSIdEnum);
207  InputUpdateFromVectorx(femmodel,N_gia_rate,NGiaRateEnum,VertexSIdEnum);
208  InputUpdateFromVectorx(femmodel,RSLg_rate,SealevelRSLRateEnum,VertexSIdEnum);
209  InputUpdateFromVectorx(femmodel,U_esa,SealevelUEsaEnum,VertexSIdEnum);
210  InputUpdateFromVectorx(femmodel,N_esa,SealevelNEsaEnum,VertexSIdEnum);
211  InputUpdateFromVectorx(femmodel,U_gia,UGiaEnum,VertexSIdEnum);
212  InputUpdateFromVectorx(femmodel,N_gia,NGiaEnum,VertexSIdEnum);
213  InputUpdateFromVectorx(femmodel,RSLg,SealevelRSLEnum,VertexSIdEnum);
214  InputUpdateFromVectorx(femmodel,RSLg_eustatic,SealevelRSLEustaticEnum,VertexSIdEnum);
215 
216  if (horiz){
217  InputUpdateFromVectorx(femmodel,U_north_esa,SealevelUNorthEsaEnum,VertexSIdEnum); // north motion
218  InputUpdateFromVectorx(femmodel,U_east_esa,SealevelUEastEsaEnum,VertexSIdEnum); // east motion
219  } /*}}}*/
220  }
221 
222  if(iscoupler){
223  /*transfer sea level back to ice caps:*/
224  TransferSealevel(femmodel,SealevelNEsaRateEnum);
225  TransferSealevel(femmodel,NGiaRateEnum);
226  TransferSealevel(femmodel,SealevelUEsaRateEnum);
227  TransferSealevel(femmodel,UGiaRateEnum);
228 
229  //reset cumdeltathickness to 0:
230  InputUpdateFromConstantx(femmodel->inputs2,femmodel->elements,0.,SealevelriseCumDeltathicknessEnum);
231  }
232 
233  /*reset counter to 1:*/
234  femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter.
235 
236  /*free ressources:{{{*/
237  delete RSLg;
238  delete RSLg_rate;
239  delete RSLg_eustatic;
240  delete U_esa;
241  delete U_esa_rate;
242  delete N_esa;
243  delete N_esa_rate;
244  delete BPg;
245 
246  if(horiz){
247  delete U_north_esa;
248  delete U_east_esa;
249  }
250  delete N_gia;
251  delete U_gia;
252  delete N_gia_rate;
253  delete U_gia_rate;
254  //delete masks;
255  /*}}}*/
256 
257 }
258 /*}}}*/
259 void dynstr_core(FemModel* femmodel){ /*{{{*/
260 
261  /*variables:*/
262  Vector<IssmDouble> *bedrock = NULL;
263  Vector<IssmDouble> *SL = NULL;
264  Vector<IssmDouble> *steric_rate_g = NULL;
265  Vector<IssmDouble> *dynamic_rate_g = NULL;
266  Vector<IssmDouble> *U_esa_rate= NULL;
267  Vector<IssmDouble> *N_esa_rate= NULL;
268  Vector<IssmDouble> *U_gia_rate= NULL;
269  Vector<IssmDouble> *N_gia_rate= NULL;
270 
271  /*parameters: */
272  bool isslr=0;
273  int solution_type;
274  IssmDouble dt;
275  int geodetic=0;
276 
277  /*Retrieve parameters:*/
278  femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
279  femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
280  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
281  femmodel->parameters->FindParam(&geodetic,SolidearthSettingsComputesealevelchangeEnum);
282 
283  /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
284  if(solution_type==SealevelriseSolutionEnum)isslr=1;
285 
286  /*Should we be here?:*/
287  if(!isslr)return;
288 
289  /*Verbose: */
290  if(VerboseSolution()) _printf0_(" computing steric sea level rise\n");
291 
292  /*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/
293  GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum);
294  GetVectorFromInputsx(&SL,femmodel,SealevelEnum,VertexSIdEnum);
295  GetStericRate(&steric_rate_g,femmodel);
296  GetDynamicRate(&dynamic_rate_g,femmodel);
297  if(geodetic){
298  GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum);
299  GetVectorFromInputsx(&U_gia_rate,femmodel,UGiaRateEnum,VertexSIdEnum);
300  GetVectorFromInputsx(&N_esa_rate,femmodel,SealevelNEsaRateEnum,VertexSIdEnum);
301  GetVectorFromInputsx(&N_gia_rate,femmodel,NGiaRateEnum,VertexSIdEnum);
302  }
303 
304  /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate) * dt + steric_rate + dynamic_rate dt*/
305  if(geodetic){
306  SL->AXPY(N_gia_rate,dt);
307  SL->AXPY(N_esa_rate,dt);
308  }
309  SL->AXPY(steric_rate_g,dt);
310  SL->AXPY(dynamic_rate_g,dt);
311 
312  /*compute new bedrock position: */
313  if(geodetic){
314  bedrock->AXPY(U_esa_rate,dt);
315  bedrock->AXPY(U_gia_rate,dt);
316  }
317 
318  /*update element inputs:*/
319  InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum);
320  InputUpdateFromVectorx(femmodel,SL,SealevelEnum,VertexSIdEnum);
321 
322  /*Free ressources:*/
323  delete bedrock;
324  delete SL;
325  delete steric_rate_g;
326  delete dynamic_rate_g;
327  if(geodetic){
328  delete U_esa_rate;
329  delete U_gia_rate;
330  delete N_esa_rate;
331  delete N_gia_rate;
332  }
333 }
334 /*}}}*/
335 
336 SealevelMasks* sealevelrise_core_masks(FemModel* femmodel) { /*{{{*/
337 
338  if(VerboseSolution()) _printf0_(" computing sea level masks\n");
339 
340  /*initialize SealevelMasks structure: */
341  SealevelMasks* masks=new SealevelMasks(femmodel->elements->Size());
342 
343  /*go through elements and fill the masks: */
344  for (int i=0;i<femmodel->elements->Size();i++){
345  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
346  element->SetSealevelMasks(masks);
347  }
348 
349  return masks;
350 }/*}}}*/
351 void sealevelrise_core_geometry(FemModel* femmodel) { /*{{{*/
352 
353  /*Geometry core where we compute indices into tables pre computed in the SealevelRiseAnalysis: */
354 
355  /*parameters: */
356  bool spherical=true;
357  IssmDouble *latitude = NULL;
358  IssmDouble *longitude = NULL;
359  IssmDouble *radius = NULL;
360  IssmDouble *xx = NULL;
361  IssmDouble *yy = NULL;
362  IssmDouble *zz = NULL;
363  int horiz;
364  bool geometrydone = false;
365 
366 
367  /*retrieve parameters:*/
368  femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
369  femmodel->parameters->FindParam(&geometrydone,SealevelriseGeometryDoneEnum);
370 
371  if(geometrydone){
372  if(VerboseSolution()) _printf0_(" geometrical offsets have already been computed, skipping \n");
373  return; //don't need to run this again.
374  }
375 
376  /*Verbose: */
377  if(VerboseSolution()) _printf0_(" computing geometrical offsets into precomputed Green tables \n");
378 
379  /*first, recover lat,long and radius vectors from vertices: */
380  VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
381  if(horiz) VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
382 
383 
384  /*Run sealevelrie geometry routine in elements:*/
385  for(int i=0;i<femmodel->elements->Size();i++){
386  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
387  element->SealevelriseGeometry(latitude,longitude,radius,xx,yy,zz);
388  }
389 
390  /*Free ressources:*/
391  if(horiz){
392  xDelete<IssmDouble>(xx);
393  xDelete<IssmDouble>(yy);
394  xDelete<IssmDouble>(zz);
395  }
396  xDelete<IssmDouble>(longitude);
397  xDelete<IssmDouble>(radius);
398 
399  /*Record the fact that we ran this module already: */
400  femmodel->parameters->SetParam(true,SealevelriseGeometryDoneEnum);
401 
402 
403 }/*}}}*/
404 Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* poceanarea){ /*{{{*/
405 
406  /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/
407 
408  Vector<IssmDouble> *RSLgi = NULL;
409  IssmDouble RSLgi_oceanaverage = 0;
410 
411  /*parameters: */
412  int gsize;
413  IssmDouble oceanarea;
414 
415  /*outputs:*/
416  IssmDouble eustatic;
417 
418  if(VerboseSolution()) _printf0_(" computing eustatic components on ice\n");
419 
420 
421  /*Figure out size of g-set deflection vector and allocate solution vector: */
422  gsize = femmodel->nodes->NumberOfDofs(GsetEnum);
423 
424  /*Initialize:*/
425  RSLgi = new Vector<IssmDouble>(gsize);
426 
427  /*call the eustatic main module: */
428  femmodel->SealevelriseEustatic(RSLgi,&oceanarea,&eustatic, masks); //this computes
429 
430  /*we need to average RSLgi over the ocean: RHS term 4 in Eq.4 of Farrel and clarke. Only the elements can do that: */
431  RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi,masks, oceanarea);
432 
433  /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the
434  * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/
435  RSLgi->Shift(-eustatic-RSLgi_oceanaverage);
436 
437  /*save eustatic value for results: */
438  femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelRSLEustaticEnum,-eustatic));
439 
440  /*Assign output pointers and return: */
441  *poceanarea=oceanarea;
442  return RSLgi;
443 }/*}}}*/
444 Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector<IssmDouble>* RSLg_eustatic,IssmDouble oceanarea){ /*{{{*/
445 
446  /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
447  non eustatic core of the SLR solution */
448 
449  Vector<IssmDouble> *RSLg = NULL;
450  Vector<IssmDouble> *RSLg_old = NULL;
451  Vector<IssmDouble> *BPg = NULL;
452 
453  Vector<IssmDouble> *RSLgo = NULL; //ocean convolution of the perturbation to gravity potential.
454  Vector<IssmDouble> *RSLgo_rot= NULL; // rotational feedback
455  IssmDouble RSLgo_oceanaverage = 0; //average of RSLgo over the ocean.
456 
457  /*parameters: */
458  int count;
459  bool save_results;
460  int gsize;
461  bool converged=true;
462  bool rotation=true;
463  bool verboseconvolution=true;
464  int max_nonlinear_iterations;
465  IssmDouble eps_rel;
466  IssmDouble eps_abs;
467  IssmDouble eustatic;
468  IssmDouble Ixz, Iyz, Izz;
469  int bp_compute_fingerprints= 0;
470 
471  if(VerboseSolution()) _printf0_(" converging on ocean components\n");
472 
473  /*Recover some parameters: */
474  femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum);
475  femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
476  femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
477 
478  /*computational flag: */
479  femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
480 
481  /*Figure out size of g-set deflection vector and allocate solution vector: */
482  gsize = femmodel->nodes->NumberOfDofs(GsetEnum);
483 
484  /*Initialize:*/
485  RSLg = new Vector<IssmDouble>(gsize);
486  RSLg->Assemble();
487  RSLg_eustatic->Copy(RSLg); //first initialize RSLg with the eustatic component computed in sealevelrise_core_eustatic.
488 
489  RSLg_old = new Vector<IssmDouble>(gsize);
490  RSLg_old->Assemble();
491 
492  count=1;
493  converged=false;
494 
495  /*Start loop: */
496  for(;;){
497 
498  //save pointer to old sea level rise
499  delete RSLg_old; RSLg_old=RSLg;
500 
501  /*Initialize solution vector: */
502  RSLg = new Vector<IssmDouble>(gsize); RSLg->Assemble();
503  RSLgo = new Vector<IssmDouble>(gsize); RSLgo->Assemble();
504 
505  /*call the non eustatic module: */
506  femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, masks, verboseconvolution);
507 
508  /*assemble solution vector: */
509  RSLgo->Assemble();
510 
511  if(rotation){
512 
513  /*call rotational feedback module: */
514  RSLgo_rot = new Vector<IssmDouble>(gsize); RSLgo_rot->Assemble();
515  femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, masks);
516  RSLgo_rot->Assemble();
517 
518  /*save changes in inertia tensor as results: */
519  femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz));
520  femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,Iyz));
521  femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorZZEnum,Izz));
522 
523  RSLgo->AXPY(RSLgo_rot,1);
524  }
525 
526  /*we need to average RSLgo over the ocean: RHS term 5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
527  RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo,masks, oceanarea);
528 
529  /*RSLg is the sum of the eustatic term, and the ocean terms: */
530  RSLg_eustatic->Copy(RSLg); RSLg->AXPY(RSLgo,1);
531  RSLg->Shift(-RSLgo_oceanaverage);
532 
533  /*convergence criterion:*/
534  slrconvergence(&converged,RSLg,RSLg_old,eps_rel,eps_abs);
535 
536 
537  /*free ressources: */
538  delete RSLgo;
539  delete RSLgo_rot;
540 
541  /*Increase count: */
542  count++;
543  if(converged==true){
544  break;
545  }
546  if(count>=max_nonlinear_iterations){
547  _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");
548  converged=true;
549  break;
550  }
551 
552  /*some minor verbosing adjustment:*/
553  if(count>1)verboseconvolution=false;
554 
555  }
556  if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n");
557 
558 
559  /*if we had bottom pressure loading, add dynamic sea level
560  * to RSL:*/
561  femmodel->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
562  if(bp_compute_fingerprints){
563  GetVectorFromInputsx(&BPg,femmodel,DslSeaWaterPressureChangeAtSeaFloor,VertexSIdEnum);
564  RSLg->AXPY(BPg,1);
565  }
566 
567  delete RSLg_old;
568  delete BPg;
569 
570  return RSLg;
571 } /*}}}*/
572 void sealevelrise_core_elastic(Vector<IssmDouble>** pU_esa, Vector<IssmDouble>** pU_north_esa,Vector<IssmDouble>** pU_east_esa,FemModel* femmodel,Vector<IssmDouble>* RSLg, SealevelMasks* masks){ /*{{{*/
573 
574  Vector<IssmDouble> *U_esa = NULL;
575  Vector<IssmDouble> *U_north_esa = NULL;
576  Vector<IssmDouble> *U_east_esa = NULL;
577 
578  /*parameters: */
579  int gsize;
580  bool spherical=true;
581 
582  IssmDouble *latitude = NULL;
583  IssmDouble *longitude = NULL;
584  IssmDouble *radius = NULL;
585  IssmDouble *xx = NULL;
586  IssmDouble *yy = NULL;
587  IssmDouble *zz = NULL;
588  int horiz;
589 
590  if(VerboseSolution()) _printf0_(" computing vertical and horizontal geodetic signatures\n");
591 
592  /*retrieve some parameters:*/
593  femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
594 
595  /*find size of vectors:*/
596  gsize = femmodel->nodes->NumberOfDofs(GsetEnum);
597 
598  /*intialize vectors:*/
599  U_esa = new Vector<IssmDouble>(gsize);
600  if (horiz){
601  U_north_esa = new Vector<IssmDouble>(gsize);
602  U_east_esa = new Vector<IssmDouble>(gsize);
603  }
604 
605  /*retrieve geometric information: */
606  VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
607  VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
608 
609  /*call the elastic main modlule:*/
610  femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg, masks);
611 
612  /*Assign output pointers:*/
613  *pU_esa=U_esa;
614  if(horiz){
615  *pU_east_esa=U_east_esa;
616  *pU_north_esa=U_north_esa;
617  }
618 
619  /*Free ressources: */
620  xDelete<IssmDouble>(longitude);
621  xDelete<IssmDouble>(latitude);
622  xDelete<IssmDouble>(xx);
623  xDelete<IssmDouble>(yy);
624  xDelete<IssmDouble>(zz);
625  xDelete<IssmDouble>(radius);
626 }
627 /*}}}*/
628 void sealevelrise_core_viscous(Vector<IssmDouble>** pU_gia, Vector<IssmDouble>** pN_gia,FemModel* femmodel,Vector<IssmDouble>* RSLg){ /*{{{*/
629 
630  /*variables:*/
631  Vector<IssmDouble> *U_gia = NULL;
632  Vector<IssmDouble> *N_gia = NULL;
633 
634  /*parameters:*/
635  int frequency;
636  IssmDouble dt;
637 
638  if(VerboseSolution()) _printf0_(" computing viscous components\n");
639 
640  /*retrieve some parameters:*/
641  femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
642  femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
643 
644  /*recover GIA rates:*/
645  GetVectorFromInputsx(&U_gia,femmodel,UGiaRateEnum,VertexSIdEnum);
646  GetVectorFromInputsx(&N_gia,femmodel,NGiaRateEnum,VertexSIdEnum);
647 
648  /*we just loaded rates, that's not what's being asked, scale by time:*/
649  U_gia->Scale(frequency*dt);
650  N_gia->Scale(frequency*dt);
651 
652  /*Assign output pointers:*/
653  *pU_gia=U_gia;
654  *pN_gia=N_gia;
655 
656 }
657 /*}}}*/
658 void sealevelrise_diagnostics(FemModel* femmodel,Vector<IssmDouble>* RSLg){ /*{{{*/
659 
660  /*compute spherical harmonics deg 1 and deg 2 coefficeints:*/
661 
662 }
663 /*}}}*/
664 void GetDynamicRate(Vector<IssmDouble> ** pdynamic_rate_g, FemModel* femmodel){ /*{{{*/
665 
666  int dslmodel=-1;
667  IssmDouble time;
668 
669  /*variables:*/
670  Vector<IssmDouble> *dynamic_rate_g = NULL;
671 
672  /*Update steric rates before retrieving them on Vertex SID set:*/
673  femmodel->parameters->FindParam(&dslmodel,DslModelEnum);
674  femmodel->parameters->FindParam(&time,TimeEnum);
675  if(dslmodel==1){
676  TransientInput2* transient_input = femmodel->inputs2->GetTransientInput(DslSeaSurfaceHeightChangeAboveGeoidEnum);
677  TriaInput2* tria_input=transient_input->GetTriaInput(time);
678  Input2* tria_input_copy=tria_input->copy();
679  tria_input_copy->ChangeEnum(DslDynamicRateEnum);
680  femmodel->inputs2->AddInput(tria_input_copy);
681  }
682  else if(dslmodel==2){
683 
684  int modelid;
685 
686  /*Recover modelid:*/
687  femmodel->parameters->FindParam(&modelid,DslModelidEnum);
688  modelid--; //from matlab.
689 
690  /*find the DslSeaSurfaceHeightChangeAboveGeoidEnum dataset of transient inputs:*/
691  DatasetInput2* dataset_input=femmodel->inputs2->GetDatasetInput2(DslSeaSurfaceHeightChangeAboveGeoidEnum);
692 
693  /*Go find the modelid'th transient input:*/
694  TriaInput2* tria_input=dataset_input->GetTriaInputByOffset(modelid);
695 
696  /*Plug into DslDynamicRate input: */
697  Input2* tria_input_copy=tria_input->copy();
698  tria_input_copy->ChangeEnum(DslDynamicRateEnum);
699  femmodel->inputs2->AddInput(tria_input_copy);
700  }
701  else _error_("not implemented yet");
702 
703  GetVectorFromInputsx(&dynamic_rate_g,femmodel,DslDynamicRateEnum,VertexSIdEnum);
704  *pdynamic_rate_g=dynamic_rate_g;
705 }
706 /*}}}*/
707 void GetStericRate(Vector<IssmDouble> ** psteric_rate_g, FemModel* femmodel){ /*{{{*/
708 
709  int dslmodel=-1;
710  IssmDouble time;
711 
712  /*variables:*/
713  Vector<IssmDouble> *steric_rate_g = NULL;
714 
715  /*Update steric rates before retrieving them on Vertex SID set:*/
716  femmodel->parameters->FindParam(&dslmodel,DslModelEnum);
717  femmodel->parameters->FindParam(&time,TimeEnum);
718  if(dslmodel==1){
719  TransientInput2* transient_input = femmodel->inputs2->GetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum);
720  TriaInput2* tria_input=transient_input->GetTriaInput(time);
721  Input2* tria_input_copy=tria_input->copy();
722  tria_input_copy->ChangeEnum(DslStericRateEnum);
723  femmodel->inputs2->AddInput(tria_input_copy);
724  }
725  else if (dslmodel==2){
726  int modelid;
727 
728  /*Recover modelid:*/
729  femmodel->parameters->FindParam(&modelid,DslModelidEnum);
730 
731  modelid--; //from matlab.
732 
733  /*find the DslGlobalAverageThermostericSeaLevelChangeEnum dataset of transient inputs:*/
734  DatasetInput2* dataset_input=femmodel->inputs2->GetDatasetInput2(DslGlobalAverageThermostericSeaLevelChangeEnum);
735 
736  /*Go find the modelid'th transient input:*/
737  TriaInput2* tria_input=dataset_input->GetTriaInputByOffset(modelid);
738 
739  /*Plug into DslStericRate input: */
740  Input2* tria_input_copy=tria_input->copy();
741  tria_input_copy->ChangeEnum(DslStericRateEnum);
742  femmodel->inputs2->AddInput(tria_input_copy);
743  }
744  else _error_("not implemented yet");
745 
746  GetVectorFromInputsx(&steric_rate_g,femmodel,DslStericRateEnum,VertexSIdEnum);
747  *psteric_rate_g=steric_rate_g;
748 }
749 /*}}}*/
750 
751 /*support routines:*/
752 void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/
753 
754  /*forcing being transferred from models to earth: */
755  IssmDouble** forcings=NULL;
756  IssmDouble* forcing=NULL;
757  Vector<IssmDouble>* forcingglobal=NULL;
758  int* nvs=NULL;
759 
760  /*transition vectors:*/
761  IssmDouble** transitions=NULL;
762  int ntransitions;
763  int* transitions_m=NULL;
764  int* transitions_n=NULL;
765  int nv;
766 
767  /*communicators:*/
768  ISSM_MPI_Comm tocomm;
769  ISSM_MPI_Comm* fromcomms=NULL;
770  ISSM_MPI_Status status;
771  int my_rank;
772  int modelid,earthid;
773  int nummodels;
774 
775  /*Recover some parameters: */
776  femmodel->parameters->FindParam(&modelid,ModelIdEnum);
777  femmodel->parameters->FindParam(&earthid,EarthIdEnum);
778  femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
779  my_rank=IssmComm::GetRank();
780 
781  /*retrieve the inter communicators that will be used to send data from each ice cap to the earth: */
782  if(modelid==earthid){
783  GenericParam<ISSM_MPI_Comm*>* parcoms = dynamic_cast<GenericParam<ISSM_MPI_Comm*>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
784  if(!parcoms)_error_("TransferForcing error message: could not find IcecapToEarthComm communicator");
785  fromcomms=parcoms->GetParameterValue();
786  }
787  else {
788  GenericParam<ISSM_MPI_Comm>* parcom = dynamic_cast<GenericParam<ISSM_MPI_Comm>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
789  if(!parcom)_error_("TransferForcing error message: could not find IcecapToEarthComm communicator");
790  tocomm=parcom->GetParameterValue();
791  }
792 
793  /*For each icecap, retrieve the forcing vector that will be sent to the earth model: */
794  if(modelid!=earthid){
795  nv=femmodel->vertices->NumberOfVertices();
796  GetVectorFromInputsx(&forcing,femmodel,forcingenum,VertexSIdEnum);
797  }
798 
799  /*Send the forcing to the earth model:{{{*/
800  if(my_rank==0){
801  if(modelid==earthid){
802  forcings=xNew<IssmDouble*>(nummodels-1);
803  nvs=xNew<int>(nummodels-1);
804  for(int i=0;i<earthid;i++){
805  ISSM_MPI_Recv(nvs+i, 1, ISSM_MPI_INT, 0,i, fromcomms[i], &status);
806  forcings[i]=xNew<IssmDouble>(nvs[i]);
807  ISSM_MPI_Recv(forcings[i], nvs[i], ISSM_MPI_DOUBLE, 0,i, fromcomms[i], &status);
808  }
809 
810  }
811  else{
812  ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, modelid, tocomm);
813  ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, tocomm);
814  }
815  }
816  /*}}}*/
817 
818  /*On the earth model, consolidate all the forcings into one, and update the elements dataset accordingly: {{{*/
819  if(modelid==earthid){
820 
821  /*Out of all the delta thicknesses, build one delta thickness vector made of all the ice cap contributions.
822  *First, build the global delta thickness vector in the earth model: */
823  nv=femmodel->vertices->NumberOfVertices();
824  GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum);
825 
826  /*Retrieve transition vectors, used to plug from each ice cap into the global forcing:*/
827  femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);
828 
829  if(ntransitions!=earthid)_error_("TransferForcing error message: number of transition vectors is not equal to the number of icecaps!");
830 
831  /*Go through all the delta thicknesses coming from each ice cap: */
832  if(my_rank==0){
833  for(int i=0;i<earthid;i++){
834 
835  IssmDouble* forcingfromcap= forcings[i]; //careful, this only exists on rank 0 of the earth model!
836  IssmDouble* transition=transitions[i];
837  int M=transitions_m[i];
838 
839  /*build index to plug values: */
840  int* index=xNew<int>(M); for(int i=0;i<M;i++)index[i]=reCast<int>(transition[i])-1; //matlab indexing!
841 
842  /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular
843  * ice cap: */
844  forcingglobal->SetValues(M,index,forcingfromcap,ADD_VAL);
845  xDelete<int>(index);
846  }
847  }
848 
849  /*Assemble vector:*/
850  forcingglobal->Assemble();
851 
852  /*Plug into elements:*/
853  InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum);
854  }
855  /*}}}*/
856 
857  /*Free ressources:{{{*/
858  if(forcings){
859  for(int i=0;i<nummodels-1;i++){
860  IssmDouble* temp=forcings[i]; xDelete<IssmDouble>(temp);
861  }
862  xDelete<IssmDouble*>(forcings);
863  }
864  if(forcing)xDelete<IssmDouble>(forcing);
865  if(forcingglobal)delete forcingglobal;
866  if(transitions){
867  for(int i=0;i<earthid;i++){
868  IssmDouble* temp=transitions[i];
869  xDelete<IssmDouble>(temp);
870  }
871  xDelete<IssmDouble*>(transitions);
872  xDelete<int>(transitions_m);
873  xDelete<int>(transitions_n);
874  }
875  if(nvs)xDelete<int>(nvs);
876  /*}}}*/
877 
878 } /*}}}*/
879 void TransferSealevel(FemModel* femmodel,int forcingenum){ /*{{{*/
880 
881  /*forcing being transferred from earth to ice caps: */
882  IssmDouble* forcing=NULL;
883  IssmDouble* forcingglobal=NULL;
884 
885  /*transition vectors:*/
886  IssmDouble** transitions=NULL;
887  int ntransitions;
888  int* transitions_m=NULL;
889  int* transitions_n=NULL;
890  int nv;
891 
892  /*communicators:*/
893  ISSM_MPI_Comm fromcomm;
894  ISSM_MPI_Comm* tocomms=NULL;
895  ISSM_MPI_Status status;
896  int my_rank;
897  int modelid,earthid;
898  int nummodels;
899  int numcoms;
900 
901  /*Recover some parameters: */
902  femmodel->parameters->FindParam(&modelid,ModelIdEnum);
903  femmodel->parameters->FindParam(&earthid,EarthIdEnum);
904  femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
905  my_rank=IssmComm::GetRank();
906 
907  /*retrieve the inter communicators that will be used to send data from earth to ice caps:*/
908  if(modelid==earthid){
909  GenericParam<ISSM_MPI_Comm*>* parcoms = dynamic_cast<GenericParam<ISSM_MPI_Comm*>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
910  if(!parcoms)_error_("TransferSealevel error message: could not find IcecapToEarthComm communicator");
911  tocomms=parcoms->GetParameterValue();
912  //femmodel->parameters->FindParam((int**)(&tocomms),&numcoms,IcecapToEarthCommEnum);
913  }
914  else{
915  GenericParam<ISSM_MPI_Comm>* parcom = dynamic_cast<GenericParam<ISSM_MPI_Comm>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
916  if(!parcom)_error_("TransferSealevel error message: could not find IcecapToEarthComm communicator");
917  fromcomm=parcom->GetParameterValue();
918  //femmodel->parameters->FindParam((int*)(&fromcomm), IcecapToEarthCommEnum);
919  }
920 
921  /*Retrieve sea-level on earth model: */
922  if(modelid==earthid){
923  nv=femmodel->vertices->NumberOfVertices();
924  GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum);
925  }
926 
927  /*Send the forcing to the ice caps:{{{*/
928  if(my_rank==0){
929 
930  if(modelid==earthid){
931 
932  /*Retrieve transition vectors, used to figure out global forcing contribution to each ice cap's own elements: */
933  femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);
934 
935  if(ntransitions!=earthid)_error_("TransferSeaLevel error message: number of transition vectors is not equal to the number of icecaps!");
936 
937  for(int i=0;i<earthid;i++){
938  nv=transitions_m[i];
939  forcing=xNew<IssmDouble>(nv);
940  IssmDouble* transition=transitions[i];
941  for(int j=0;j<nv;j++){
942  forcing[j]=forcingglobal[reCast<int>(transition[j])-1];
943  }
944  ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, i, tocomms[i]);
945  ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, i, tocomms[i]);
946  }
947  }
948  else{
949  ISSM_MPI_Recv(&nv, 1, ISSM_MPI_INT, 0, modelid, fromcomm, &status);
950  forcing=xNew<IssmDouble>(nv);
951  ISSM_MPI_Recv(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, fromcomm, &status);
952  }
953  }
954  /*}}}*/
955 
956  /*On each ice cap, spread the forcing across cpus, and update the elements dataset accordingly: {{{*/
957  if(modelid!=earthid){
958 
959  ISSM_MPI_Bcast(&nv,1,ISSM_MPI_INT,0,IssmComm::GetComm());
960  if(my_rank!=0)forcing=xNew<IssmDouble>(nv);
961  ISSM_MPI_Bcast(forcing,nv,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
962 
963  /*Plug into elements:*/
964  InputUpdateFromVectorx(femmodel,forcing,forcingenum,VertexSIdEnum);
965  }
966  /*}}}*/
967 
968  /*Free ressources:{{{*/
969  if(forcingglobal)xDelete<IssmDouble>(forcingglobal);
970  if(forcing)xDelete<IssmDouble>(forcing);
971  if(transitions){
972  for(int i=0;i<ntransitions;i++){
973  IssmDouble* temp=transitions[i];
974  xDelete<IssmDouble>(temp);
975  }
976  xDelete<IssmDouble*>(transitions);
977  xDelete<int>(transitions_m);
978  xDelete<int>(transitions_n);
979  }
980  /*}}}*/
981 
982 } /*}}}*/
983 void EarthMassTransport(FemModel* femmodel){ /*{{{*/
984 
985  IssmDouble time,dt;
986  Vector<IssmDouble> *oldthickness = NULL;
987  Vector<IssmDouble> *newthickness = NULL;
988  Vector<IssmDouble> *deltathickness = NULL;
989  Vector<IssmDouble> *cumdeltathickness = NULL;
990  int nv;
991 
992  if(VerboseSolution()) _printf0_(" computing earth mass transport\n");
993 
994  /*This module has to be recoded! We do not grab spc thicknesses from an earth model
995  * anymore! Thicknesses come from a mass transport module applied to each basin!
996  * Commeting out for now:*/
997  _error_("EarthMassTransport error message: not supported anymore!");
998 
999  /*This mass transport module for the Earth is because we might have thickness variations as spcs
1000  * specified in the md.slr class, outside of what we will get from the icecaps. That's why we get t
1001  * the thickness variations from SealevelriseSpcthicknessEnum.*/
1002 
1003  /*No mass transport module was called, so we are just going to retrieve the geometry thickness
1004  * at this time step, at prior time step, and plug the difference as deltathickness: */
1005  /*femmodel->parameters->FindParam(&time,TimeEnum);
1006  femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
1007  nv=femmodel->vertices->NumberOfVertices();
1008 
1009  GetVectorFromInputsx(&newthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum);
1010  GetVectorFromInputsx(&oldthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum,time-dt);*/
1011 
1012  /*compute deltathickness: */
1013  /*deltathickness = new Vector<IssmDouble>(nv);
1014  newthickness->Copy(deltathickness); deltathickness->AXPY(oldthickness,-1); */
1015 
1016  /*plug into elements:*/
1017 // InputUpdateFromVectorx(femmodel,deltathickness,SurfaceloadIceThicknessChangeEnum,VertexSIdEnum);
1018 
1019  /*add to cumulated delta thickness: */
1020  /*GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
1021  cumdeltathickness->AXPY(deltathickness,1);
1022  InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);*/
1023 
1024  /*free ressources:*/
1025  /*delete oldthickness;
1026  delete newthickness;
1027  delete deltathickness;
1028  delete cumdeltathickness;*/
1029 
1030 } /*}}}*/
1031 void slrconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/
1032 
1033  bool converged=true;
1034  IssmDouble ndS,nS;
1035  Vector<IssmDouble> *dRSLg = NULL;
1036 
1037  //compute norm(du) and norm(u) if requested
1038  dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0);
1039  ndS=dRSLg->Norm(NORM_TWO);
1040 
1041  if (xIsNan<IssmDouble>(ndS)) _error_("convergence criterion is NaN!");
1042 
1043  if(!xIsNan<IssmDouble>(eps_rel)){
1044  nS=RSLg_old->Norm(NORM_TWO);
1045  if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN!");
1046  }
1047 
1048  //clean up
1049  delete dRSLg;
1050 
1051  //print
1052  if(!xIsNan<IssmDouble>(eps_rel)){
1053  if((ndS/nS)<eps_rel){
1054  if(VerboseConvergence()) _printf0_(setw(50) << left << " convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " < " << eps_rel*100 << " %\n");
1055  }
1056  else{
1057  if(VerboseConvergence()) _printf0_(setw(50) << left << " convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " > " << eps_rel*100 << " %\n");
1058  converged=false;
1059  }
1060  }
1061  if(!xIsNan<IssmDouble>(eps_abs)){
1062  if(ndS<eps_abs){
1063  if(VerboseConvergence()) _printf0_(setw(50) << left << " convergence criterion: norm(dS)" << ndS << " < " << eps_abs << " \n");
1064  }
1065  else{
1066  if(VerboseConvergence()) _printf0_(setw(50) << left << " convergence criterion: norm(dS)" << ndS << " > " << eps_abs << " \n");
1067  converged=false;
1068  }
1069  }
1070 
1071  /*assign output*/
1072  *pconverged=converged;
1073 
1074 } /*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
DslSeaWaterPressureChangeAtSeaFloor
@ DslSeaWaterPressureChangeAtSeaFloor
Definition: EnumDefinitions.h:542
DslSeaSurfaceHeightChangeAboveGeoidEnum
@ DslSeaSurfaceHeightChangeAboveGeoidEnum
Definition: EnumDefinitions.h:541
DslGlobalAverageThermostericSeaLevelChangeEnum
@ DslGlobalAverageThermostericSeaLevelChangeEnum
Definition: EnumDefinitions.h:540
SaveResultsEnum
@ SaveResultsEnum
Definition: EnumDefinitions.h:302
SealevelUEsaRateEnum
@ SealevelUEsaRateEnum
Definition: EnumDefinitions.h:686
SolidearthSettingsAbstolEnum
@ SolidearthSettingsAbstolEnum
Definition: EnumDefinitions.h:306
SealevelUEastEsaEnum
@ SealevelUEastEsaEnum
Definition: EnumDefinitions.h:684
FemModel::vertices
Vertices * vertices
Definition: FemModel.h:49
TransientIsmasstransportEnum
@ TransientIsmasstransportEnum
Definition: EnumDefinitions.h:449
IssmDouble
double IssmDouble
Definition: types.h:37
DatasetInput2
Definition: DatasetInput2.h:14
sealevelrise_core_masks
SealevelMasks * sealevelrise_core_masks(FemModel *femmodel)
Definition: sealevelchange_core.cpp:336
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
DslModelidEnum
@ DslModelidEnum
Definition: EnumDefinitions.h:126
VerboseConvergence
bool VerboseConvergence(void)
Definition: Verbosity.cpp:26
SealevelUEsaEnum
@ SealevelUEsaEnum
Definition: EnumDefinitions.h:685
SealevelUNorthEsaEnum
@ SealevelUNorthEsaEnum
Definition: EnumDefinitions.h:687
cores.h
FemModel::nummodels
int nummodels
Definition: FemModel.h:39
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
TimeEnum
@ TimeEnum
Definition: EnumDefinitions.h:427
ADD_VAL
@ ADD_VAL
Definition: toolkitsenums.h:14
sealevelrise_diagnostics
void sealevelrise_diagnostics(FemModel *femmodel, Vector< IssmDouble > *RSLg)
Definition: sealevelchange_core.cpp:658
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
BedEnum
@ BedEnum
Definition: EnumDefinitions.h:499
TransientInput2
Definition: TransientInput2.h:13
GetVectorFromInputsx
void GetVectorFromInputsx(IssmDouble **pvector, int *pvector_size, FemModel *femmodel, int name)
Definition: GetVectorFromInputsx.cpp:81
FemModel::results
Results * results
Definition: FemModel.h:48
Vertices::NumberOfVertices
int NumberOfVertices(void)
Definition: Vertices.cpp:255
TransientIsslrEnum
@ TransientIsslrEnum
Definition: EnumDefinitions.h:452
UGiaRateEnum
@ UGiaRateEnum
Definition: EnumDefinitions.h:594
SealevelRSLRateEnum
@ SealevelRSLRateEnum
Definition: EnumDefinitions.h:683
Vector::Duplicate
Vector< doubletype > * Duplicate(void)
Definition: Vector.h:230
SealevelNEsaEnum
@ SealevelNEsaEnum
Definition: EnumDefinitions.h:678
DslStericRateEnum
@ DslStericRateEnum
Definition: EnumDefinitions.h:543
IcecapToEarthCommEnum
@ IcecapToEarthCommEnum
Definition: EnumDefinitions.h:200
SealevelriseSolutionEnum
@ SealevelriseSolutionEnum
Definition: EnumDefinitions.h:1267
SolidearthSettingsHorizEnum
@ SolidearthSettingsHorizEnum
Definition: EnumDefinitions.h:323
UGiaEnum
@ UGiaEnum
Definition: EnumDefinitions.h:593
SealevelriseRunCountEnum
@ SealevelriseRunCountEnum
Definition: EnumDefinitions.h:331
EarthMassTransport
void EarthMassTransport(FemModel *femmodel)
Definition: sealevelchange_core.cpp:983
grd_core
void grd_core(FemModel *femmodel)
Definition: sealevelchange_core.cpp:84
Input2::ChangeEnum
void ChangeEnum(int newenumtype)
Definition: Input2.h:26
NORM_TWO
@ NORM_TWO
Definition: toolkitsenums.h:15
FemModel::RequestedDependentsx
void RequestedDependentsx(void)
Definition: FemModel.cpp:2220
Vector::Norm
doubletype Norm(NormMode norm_type)
Definition: Vector.h:342
SealevelriseTransitionsEnum
@ SealevelriseTransitionsEnum
Definition: EnumDefinitions.h:332
sealevelrise_core_geometry
void sealevelrise_core_geometry(FemModel *femmodel)
Definition: sealevelchange_core.cpp:351
Element
Definition: Element.h:41
Vector::Scale
void Scale(doubletype scale_factor)
Definition: Vector.h:355
TriaInput2::copy
Input2 * copy()
Definition: TriaInput2.cpp:72
DslModelEnum
@ DslModelEnum
Definition: EnumDefinitions.h:125
SolidearthSettingsReltolEnum
@ SolidearthSettingsReltolEnum
Definition: EnumDefinitions.h:327
SealevelriseCumDeltathicknessEnum
@ SealevelriseCumDeltathicknessEnum
Definition: EnumDefinitions.h:688
ISSM_MPI_DOUBLE
#define ISSM_MPI_DOUBLE
Definition: issmmpi.h:125
DslComputeFingerprintsEnum
@ DslComputeFingerprintsEnum
Definition: EnumDefinitions.h:128
Parameters::SetParam
void SetParam(bool boolean, int enum_type)
Definition: Parameters.cpp:441
ISSM_MPI_INT
#define ISSM_MPI_INT
Definition: issmmpi.h:127
GenericParam::GetParameterValue
P & GetParameterValue()
Definition: GenericParam.h:62
ModelIdEnum
@ ModelIdEnum
Definition: EnumDefinitions.h:274
SolidearthSettingsMaxiterEnum
@ SolidearthSettingsMaxiterEnum
Definition: EnumDefinitions.h:324
Vector::AXPY
void AXPY(Vector *X, doubletype a)
Definition: Vector.h:256
VertexSIdEnum
@ VertexSIdEnum
Definition: EnumDefinitions.h:1325
Vector::SetValues
void SetValues(int ssize, int *list, doubletype *values, InsMode mode)
Definition: Vector.h:153
FemModel::nodes
Nodes * nodes
Definition: FemModel.h:56
SealevelInertiaTensorXZEnum
@ SealevelInertiaTensorXZEnum
Definition: EnumDefinitions.h:1261
InputUpdateFromVectorx
void InputUpdateFromVectorx(FemModel *femmodel, Vector< IssmDouble > *vector, int name, int type)
Definition: InputUpdateFromVectorx.cpp:9
EarthIdEnum
@ EarthIdEnum
Definition: EnumDefinitions.h:129
TriaInput2
Definition: TriaInput2.h:8
SealevelRSLEustaticEnum
@ SealevelRSLEustaticEnum
Definition: EnumDefinitions.h:681
NGiaRateEnum
@ NGiaRateEnum
Definition: EnumDefinitions.h:592
slrconvergence
void slrconvergence(bool *pconverged, Vector< IssmDouble > *RSLg, Vector< IssmDouble > *RSLg_old, IssmDouble eps_rel, IssmDouble eps_abs)
Definition: sealevelchange_core.cpp:1031
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
sealevelrise_core_eustatic
Vector< IssmDouble > * sealevelrise_core_eustatic(FemModel *femmodel, SealevelMasks *masks, IssmDouble *poceanarea)
Definition: sealevelchange_core.cpp:404
gia_core
void gia_core(FemModel *femmodel)
Definition: gia_core.cpp:13
ISSM_MPI_Send
int ISSM_MPI_Send(void *buf, int count, ISSM_MPI_Datatype datatype, int dest, int tag, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:484
FemModel::inputs2
Inputs2 * inputs2
Definition: FemModel.h:47
sealevelrise_core_viscous
void sealevelrise_core_viscous(Vector< IssmDouble > **pU_gia, Vector< IssmDouble > **pN_gia, FemModel *femmodel, Vector< IssmDouble > *RSLg)
Definition: sealevelchange_core.cpp:628
InputDuplicatex
void InputDuplicatex(FemModel *femmodel, int original_enum, int new_enum)
Definition: InputDuplicatex.cpp:10
TransferForcing
void TransferForcing(FemModel *femmodel, int forcingenum)
Definition: sealevelchange_core.cpp:752
SolutionTypeEnum
@ SolutionTypeEnum
Definition: EnumDefinitions.h:398
ISSM_MPI_Status
int ISSM_MPI_Status
Definition: issmmpi.h:121
SealevelriseAnalysisEnum
@ SealevelriseAnalysisEnum
Definition: EnumDefinitions.h:1266
Results::AddResult
int AddResult(ExternalResult *result)
Definition: Results.cpp:33
dynstr_core
void dynstr_core(FemModel *femmodel)
Definition: sealevelchange_core.cpp:259
InputUpdateFromConstantx
void InputUpdateFromConstantx(FemModel *femmodel, bool constant, int name)
Definition: InputUpdateFromConstantx.cpp:10
TransientIscouplerEnum
@ TransientIscouplerEnum
Definition: EnumDefinitions.h:443
Profiler::Stop
void Stop(int tagenum, bool dontmpisync=true)
Definition: Profiler.cpp:179
SealevelriseGeometryDoneEnum
@ SealevelriseGeometryDoneEnum
Definition: EnumDefinitions.h:309
FemModel::solution_type
int solution_type
Definition: FemModel.h:40
DslDynamicRateEnum
@ DslDynamicRateEnum
Definition: EnumDefinitions.h:544
FemModel::elements
Elements * elements
Definition: FemModel.h:44
ISSM_MPI_Bcast
int ISSM_MPI_Bcast(void *buffer, int count, ISSM_MPI_Datatype datatype, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:162
SLRCORE
#define SLRCORE
Definition: Profiler.h:28
Vector::Assemble
void Assemble(void)
Definition: Vector.h:142
Input2
Definition: Input2.h:18
SealevelRSLEnum
@ SealevelRSLEnum
Definition: EnumDefinitions.h:680
FemModel
Definition: FemModel.h:31
SealevelEnum
@ SealevelEnum
Definition: EnumDefinitions.h:675
Inputs2::GetTransientInput
TransientInput2 * GetTransientInput(int enum_type)
Definition: Inputs2.cpp:406
GenericExternalResult
Definition: GenericExternalResult.h:21
IssmComm::GetRank
static int GetRank(void)
Definition: IssmComm.cpp:34
sealevelrise_core_elastic
void sealevelrise_core_elastic(Vector< IssmDouble > **pU_esa, Vector< IssmDouble > **pU_north_esa, Vector< IssmDouble > **pU_east_esa, FemModel *femmodel, Vector< IssmDouble > *RSLg, SealevelMasks *masks)
Definition: sealevelchange_core.cpp:572
SealevelInertiaTensorZZEnum
@ SealevelInertiaTensorZZEnum
Definition: EnumDefinitions.h:1263
Vector::AYPX
void AYPX(Vector *X, doubletype a)
Definition: Vector.h:267
sealevelchange_core
void sealevelchange_core(FemModel *femmodel)
Definition: sealevelchange_core.cpp:17
SealevelriseRequestedOutputsEnum
@ SealevelriseRequestedOutputsEnum
Definition: EnumDefinitions.h:328
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Profiler::Start
void Start(int tagenum, bool dontmpisync=true)
Definition: Profiler.cpp:139
SolidearthSettingsRotationEnum
@ SolidearthSettingsRotationEnum
Definition: EnumDefinitions.h:330
VerboseSolution
bool VerboseSolution(void)
Definition: Verbosity.cpp:24
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
SurfaceloadIceThicknessChangeEnum
@ SurfaceloadIceThicknessChangeEnum
Definition: EnumDefinitions.h:691
ISSM_MPI_Comm
int ISSM_MPI_Comm
Definition: issmmpi.h:118
SolidearthSettingsRunFrequencyEnum
@ SolidearthSettingsRunFrequencyEnum
Definition: EnumDefinitions.h:321
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
Inputs2::GetDatasetInput2
DatasetInput2 * GetDatasetInput2(int enum_type)
Definition: Inputs2.cpp:438
VertexCoordinatesx
void VertexCoordinatesx(IssmDouble **px, IssmDouble **py, IssmDouble **pz, Vertices *vertices, bool spherical)
Definition: VertexCoordinatesx.cpp:11
ISSM_MPI_Recv
int ISSM_MPI_Recv(void *buf, int count, ISSM_MPI_Datatype datatype, int source, int tag, ISSM_MPI_Comm comm, ISSM_MPI_Status *status)
Definition: issmmpi.cpp:342
Vector::Copy
void Copy(Vector *to)
Definition: Vector.h:319
TransientInput2::GetTriaInput
TriaInput2 * GetTriaInput()
Definition: TransientInput2.cpp:270
Parameters::FindParamObject
Param * FindParamObject(int enum_type)
Definition: Parameters.cpp:588
Vector::Shift
void Shift(doubletype shift)
Definition: Vector.h:309
TransferSealevel
void TransferSealevel(FemModel *femmodel, int forcingenum)
Definition: sealevelchange_core.cpp:879
Nodes::NumberOfDofs
int NumberOfDofs(int setenum)
Definition: Nodes.cpp:314
SealevelInertiaTensorYZEnum
@ SealevelInertiaTensorYZEnum
Definition: EnumDefinitions.h:1262
FemModel::profiler
Profiler * profiler
Definition: FemModel.h:42
Inputs2::AddInput
void AddInput(Input2 *in_input)
Definition: Inputs2.cpp:147
Vector< IssmDouble >
NGiaEnum
@ NGiaEnum
Definition: EnumDefinitions.h:591
GetStericRate
void GetStericRate(Vector< IssmDouble > **psteric_rate_g, FemModel *femmodel)
Definition: sealevelchange_core.cpp:707
SealevelNEsaRateEnum
@ SealevelNEsaRateEnum
Definition: EnumDefinitions.h:679
SolidearthSettingsComputesealevelchangeEnum
@ SolidearthSettingsComputesealevelchangeEnum
Definition: EnumDefinitions.h:320
GetDynamicRate
void GetDynamicRate(Vector< IssmDouble > **pdynamic_rate_g, FemModel *femmodel)
Definition: sealevelchange_core.cpp:664
SealevelMasks
Definition: SealevelMasks.h:10
DatasetInput2::GetTriaInputByOffset
TriaInput2 * GetTriaInputByOffset(int i)
Definition: DatasetInput2.cpp:222
NumModelsEnum
@ NumModelsEnum
Definition: EnumDefinitions.h:276
GenericParam
Definition: GenericParam.h:26
sealevelrise_core_noneustatic
Vector< IssmDouble > * sealevelrise_core_noneustatic(FemModel *femmodel, SealevelMasks *masks, Vector< IssmDouble > *RSLg_eustatic, IssmDouble oceanarea)
Definition: sealevelchange_core.cpp:444
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:18 for Ice Sheet System Model by   doxygen 1.8.19