sealevelchange_core.cpp File Reference

: core of the sea-level change solution More...

#include "./cores.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../classes/Inputs2/TriaInput2.h"
#include "../classes/Inputs2/TransientInput2.h"
#include "../classes/Inputs2/DatasetInput2.h"
#include "../shared/shared.h"
#include "../modules/modules.h"
#include "../solutionsequences/solutionsequences.h"

Go to the source code of this file.

Functions
void	sealevelchange_core (FemModel *femmodel)
void	grd_core (FemModel *femmodel)
void	dynstr_core (FemModel *femmodel)
SealevelMasks *	sealevelrise_core_masks (FemModel *femmodel)
void	sealevelrise_core_geometry (FemModel *femmodel)
Vector< IssmDouble > *	sealevelrise_core_eustatic (FemModel *femmodel, SealevelMasks *masks, IssmDouble *poceanarea)
Vector< IssmDouble > *	sealevelrise_core_noneustatic (FemModel *femmodel, SealevelMasks *masks, Vector< IssmDouble > *RSLg_eustatic, IssmDouble oceanarea)
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)
void	sealevelrise_core_viscous (Vector< IssmDouble > **pU_gia, Vector< IssmDouble > **pN_gia, FemModel *femmodel, Vector< IssmDouble > *RSLg)
void	sealevelrise_diagnostics (FemModel *femmodel, Vector< IssmDouble > *RSLg)
void	GetDynamicRate (Vector< IssmDouble > **pdynamic_rate_g, FemModel *femmodel)
void	GetStericRate (Vector< IssmDouble > **psteric_rate_g, FemModel *femmodel)
void	TransferForcing (FemModel *femmodel, int forcingenum)
void	TransferSealevel (FemModel *femmodel, int forcingenum)
void	EarthMassTransport (FemModel *femmodel)
void	slrconvergence (bool *pconverged, Vector< IssmDouble > *RSLg, Vector< IssmDouble > *RSLg_old, IssmDouble eps_rel, IssmDouble eps_abs)

Detailed Description

: core of the sea-level change solution

Definition in file sealevelchange_core.cpp.

Function Documentation

sealevelchange_core()

void sealevelchange_core

(

FemModel *

femmodel

)

Definition at line 17 of file sealevelchange_core.cpp.

                                            { /*{{{*/
 
   /*Start profiler*/
   femmodel->profiler->Start(SLRCORE);
 
   /*Parameters, variables:*/
   bool save_results;
   bool isslr=0;
   bool isgia=0;
   int solution_type;
 
   /*Retrieve parameters:*/
   femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
   femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
   femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
   /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
   if(solution_type==SealevelriseSolutionEnum){
      isslr=1;
      isgia=1;
   }
 
   /*Should we be here?:*/
   if(!isslr)return;
 
   /*Verbose: */
   if(VerboseSolution()) _printf0_("   computing sea level rise\n");
 
   /*Run gia core: */
   if(isgia){
      #ifdef _HAVE_GIA_
      gia_core(femmodel);
      #else
      _error_("ISSM was not compiled with gia capabilities. Exiting");
      #endif
   }
 
   /*set SLR configuration: */
   femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
 
   /*run geometry core: */
   sealevelrise_core_geometry(femmodel);
 
   /*Run geodetic:*/
   grd_core(femmodel);
 
   /*Run steric core for sure:*/
   dynstr_core(femmodel);
 
   /*Save results: */
   if(save_results){
      int        numoutputs        = 0;
      char     **requested_outputs = NULL;
 
      if(VerboseSolution()) _printf0_("   saving results\n");
      femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum);
      femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
      if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
   }
 
   /*requested dependents: */
   if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx();
   
   /*End profiler*/
   femmodel->profiler->Stop(SLRCORE);
}

grd_core()

void grd_core

(

FemModel *

femmodel

)

Definition at line 84 of file sealevelchange_core.cpp.

                                 { /*{{{*/
 
   /*Gravity rotation deformation core GRD: */
 
   /*variables:*/
   Vector<IssmDouble> *RSLg    = NULL;
   Vector<IssmDouble> *BPg    = NULL;
   Vector<IssmDouble> *RSLg_rate    = NULL;
   Vector<IssmDouble> *RSLg_eustatic  = NULL; 
   Vector<IssmDouble> *U_esa  = NULL; 
   Vector<IssmDouble> *U_esa_rate  = NULL; 
   Vector<IssmDouble> *N_esa  = NULL; 
   Vector<IssmDouble> *N_esa_rate  = NULL; 
   Vector<IssmDouble> *U_north_esa   = NULL; 
   Vector<IssmDouble> *U_east_esa    = NULL; 
   Vector<IssmDouble> *N_gia= NULL; 
   Vector<IssmDouble> *U_gia= NULL; 
   Vector<IssmDouble> *N_gia_rate= NULL; 
   Vector<IssmDouble> *U_gia_rate= NULL; 
   SealevelMasks* masks=NULL;
 
   /*parameters:*/
   bool iscoupler;
   int  solution_type;
   int  modelid,earthid;
   bool istransientmasstransport;
   int  frequency,count;
   int  horiz;
   int  geodetic=0;
   IssmDouble dt;
   IssmDouble oceanarea;
   int        bp_compute_fingerprints=0;
 
   /*Should we even be here?:*/
   femmodel->parameters->FindParam(&geodetic,SolidearthSettingsComputesealevelchangeEnum); if(!geodetic)return;
 
   /*Verbose: */
   if(VerboseSolution()) _printf0_("     computing geodetic sea level rise\n");
 
   /*retrieve more parameters:*/
   femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
   femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
   femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum);
   femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
   femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
 
   if(iscoupler){
      femmodel->parameters->FindParam(&modelid,ModelIdEnum);
      femmodel->parameters->FindParam(&earthid,EarthIdEnum);
      femmodel->parameters->FindParam(&istransientmasstransport,TransientIsmasstransportEnum);
   }
   else{
      /* we are here, we are not running in a coupler, so we will indeed compute SLR,
       * so make sure we are identified as being the Earth.:*/
      modelid=1; earthid=1; 
      /* in addition, if we are running solution_type SealevelriseSolutionEnum, make sure we 
       * run, irresepective of the time settings:*/
      count=frequency;
   }
 
   /*If we are running in coupled mode, the Earth model needs to run its own mass transport (if 
    * not already done by the mass trasnport module. For ice caps, they rely on the transient mass 
    * transport module exclusively:*/
   if(iscoupler) if(modelid==earthid) if(!istransientmasstransport) EarthMassTransport(femmodel);
 
   /*increment counter, or call solution core if count==frequency:*/
   if (count<frequency){
      count++; femmodel->parameters->SetParam(count,SealevelriseRunCountEnum); 
      return;
   }
 
   /*call sea-level rise sub cores:*/
   if(iscoupler){
      /*transfer cumulated deltathickness forcing from ice caps to earth model: */
      TransferForcing(femmodel,SealevelriseCumDeltathicknessEnum);
 
      /*we have accumulated thicknesses, dump them in deltathcikness: */
      if(modelid==earthid)InputDuplicatex(femmodel,SealevelriseCumDeltathicknessEnum,SurfaceloadIceThicknessChangeEnum);
   }
 
   /*run cores:*/
   if(modelid==earthid){
 
      /*call masks core: */
      masks=sealevelrise_core_masks(femmodel);
 
      /*call eustatic core  (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */
      RSLg_eustatic=sealevelrise_core_eustatic(femmodel,masks,&oceanarea); 
 
      /*call non-eustatic core (ocean loading tems  - 2nd and 5th terms on the RHS of Farrel and Clark) */
      RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,oceanarea); 
 
      /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */
      sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,masks);
 
      /*compute viscosus (GIA) geodetic signatures:*/
      sealevelrise_core_viscous(&U_gia,&N_gia,femmodel,RSLg);
 
      /*compute sea-level rise (low-order spherical harmonics coefficients) diagnostics:*/
      sealevelrise_diagnostics(femmodel,RSLg);
 
      /*recover N_esa  = U_esa + RSLg:*/
      N_esa=U_esa->Duplicate(); U_esa->Copy(N_esa); N_esa->AXPY(RSLg,1);
 
      /*if we had bottom pressure loading, remove dynamic sea level from geoid: */
      femmodel->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
      if(bp_compute_fingerprints){
         GetVectorFromInputsx(&BPg,femmodel,DslSeaWaterPressureChangeAtSeaFloor,VertexSIdEnum);
         N_esa->AXPY(BPg,-1);
      }
 
      /*transform these values into rates (as we only run this once each frequency turn:*/
      N_esa_rate=N_esa->Duplicate(); N_esa->Copy(N_esa_rate); N_esa_rate->Scale(1/(dt*frequency));
      U_esa_rate=U_esa->Duplicate(); U_esa->Copy(U_esa_rate); U_esa_rate->Scale(1/(dt*frequency));
      N_gia_rate=N_gia->Duplicate(); N_gia->Copy(N_gia_rate); N_gia_rate->Scale(1/(dt*frequency));
      U_gia_rate=U_gia->Duplicate(); U_gia->Copy(U_gia_rate); U_gia_rate->Scale(1/(dt*frequency));
      RSLg_rate=RSLg->Duplicate(); RSLg->Copy(RSLg_rate); RSLg_rate->Scale(1/(dt*frequency));
 
      /*get some results into elements:{{{*/
      InputUpdateFromVectorx(femmodel,U_esa_rate,SealevelUEsaRateEnum,VertexSIdEnum); 
      InputUpdateFromVectorx(femmodel,N_esa_rate,SealevelNEsaRateEnum,VertexSIdEnum); 
      InputUpdateFromVectorx(femmodel,U_gia_rate,UGiaRateEnum,VertexSIdEnum);
      InputUpdateFromVectorx(femmodel,N_gia_rate,NGiaRateEnum,VertexSIdEnum); 
      InputUpdateFromVectorx(femmodel,RSLg_rate,SealevelRSLRateEnum,VertexSIdEnum); 
      InputUpdateFromVectorx(femmodel,U_esa,SealevelUEsaEnum,VertexSIdEnum); 
      InputUpdateFromVectorx(femmodel,N_esa,SealevelNEsaEnum,VertexSIdEnum); 
      InputUpdateFromVectorx(femmodel,U_gia,UGiaEnum,VertexSIdEnum); 
      InputUpdateFromVectorx(femmodel,N_gia,NGiaEnum,VertexSIdEnum); 
      InputUpdateFromVectorx(femmodel,RSLg,SealevelRSLEnum,VertexSIdEnum); 
      InputUpdateFromVectorx(femmodel,RSLg_eustatic,SealevelRSLEustaticEnum,VertexSIdEnum); 
 
      if (horiz){
         InputUpdateFromVectorx(femmodel,U_north_esa,SealevelUNorthEsaEnum,VertexSIdEnum);   // north motion 
         InputUpdateFromVectorx(femmodel,U_east_esa,SealevelUEastEsaEnum,VertexSIdEnum);     // east motion 
      } /*}}}*/
   }
 
   if(iscoupler){
      /*transfer sea level back to ice caps:*/
      TransferSealevel(femmodel,SealevelNEsaRateEnum);
      TransferSealevel(femmodel,NGiaRateEnum);
      TransferSealevel(femmodel,SealevelUEsaRateEnum);
      TransferSealevel(femmodel,UGiaRateEnum);
 
      //reset cumdeltathickness  to 0: 
      InputUpdateFromConstantx(femmodel->inputs2,femmodel->elements,0.,SealevelriseCumDeltathicknessEnum);
   }
 
   /*reset counter to 1:*/
   femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter.
 
   /*free ressources:{{{*/
   delete RSLg;
   delete RSLg_rate;
   delete RSLg_eustatic;
   delete U_esa;
   delete U_esa_rate;
   delete N_esa;
   delete N_esa_rate;
   delete BPg;
 
   if(horiz){
      delete U_north_esa;
      delete U_east_esa;
   }
   delete N_gia;
   delete U_gia;
   delete N_gia_rate;
   delete U_gia_rate;
   //delete masks;
   /*}}}*/
 
} 

dynstr_core()

void dynstr_core

(

FemModel *

femmodel

)

Definition at line 259 of file sealevelchange_core.cpp.

                                    { /*{{{*/
 
   /*variables:*/
   Vector<IssmDouble> *bedrock  = NULL; 
   Vector<IssmDouble> *SL  = NULL; 
   Vector<IssmDouble> *steric_rate_g  = NULL; 
   Vector<IssmDouble> *dynamic_rate_g = NULL;
   Vector<IssmDouble> *U_esa_rate= NULL;
   Vector<IssmDouble> *N_esa_rate= NULL;
   Vector<IssmDouble> *U_gia_rate= NULL;
   Vector<IssmDouble> *N_gia_rate= NULL;
 
   /*parameters: */
   bool isslr=0;
   int  solution_type;
   IssmDouble          dt;
   int  geodetic=0;
 
   /*Retrieve parameters:*/
   femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
   femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
   femmodel->parameters->FindParam(&geodetic,SolidearthSettingsComputesealevelchangeEnum); 
 
   /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
   if(solution_type==SealevelriseSolutionEnum)isslr=1;
 
   /*Should we be here?:*/
   if(!isslr)return;
 
   /*Verbose: */
   if(VerboseSolution()) _printf0_("     computing steric sea level rise\n");
 
   /*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/
   GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum);
   GetVectorFromInputsx(&SL,femmodel,SealevelEnum,VertexSIdEnum);
   GetStericRate(&steric_rate_g,femmodel);
   GetDynamicRate(&dynamic_rate_g,femmodel);
   if(geodetic){
      GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum);
      GetVectorFromInputsx(&U_gia_rate,femmodel,UGiaRateEnum,VertexSIdEnum);
      GetVectorFromInputsx(&N_esa_rate,femmodel,SealevelNEsaRateEnum,VertexSIdEnum);
      GetVectorFromInputsx(&N_gia_rate,femmodel,NGiaRateEnum,VertexSIdEnum);
   }
 
   /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate)  * dt + steric_rate + dynamic_rate dt*/
   if(geodetic){
      SL->AXPY(N_gia_rate,dt);
      SL->AXPY(N_esa_rate,dt);
   }
   SL->AXPY(steric_rate_g,dt);
   SL->AXPY(dynamic_rate_g,dt);
 
   /*compute new bedrock position: */
   if(geodetic){
      bedrock->AXPY(U_esa_rate,dt);
      bedrock->AXPY(U_gia_rate,dt);
   }
 
   /*update element inputs:*/
   InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum);   
   InputUpdateFromVectorx(femmodel,SL,SealevelEnum,VertexSIdEnum);   
 
   /*Free ressources:*/ 
   delete bedrock;
   delete SL;
   delete steric_rate_g;
   delete dynamic_rate_g;
   if(geodetic){
      delete U_esa_rate;
      delete U_gia_rate;
      delete N_esa_rate;
      delete N_gia_rate;
   }
}

sealevelrise_core_masks()

SealevelMasks* sealevelrise_core_masks

(

FemModel *

femmodel

)

Definition at line 336 of file sealevelchange_core.cpp.

                                                           {  /*{{{*/
 
   if(VerboseSolution()) _printf0_("     computing sea level masks\n");
   
   /*initialize SealevelMasks structure: */
   SealevelMasks* masks=new SealevelMasks(femmodel->elements->Size());
 
   /*go through elements and fill the masks: */
   for (int i=0;i<femmodel->elements->Size();i++){
      Element*   element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      element->SetSealevelMasks(masks);
   }
 
   return masks;
}/*}}}*/

sealevelrise_core_geometry()

void sealevelrise_core_geometry

(

FemModel *

femmodel

)

Definition at line 351 of file sealevelchange_core.cpp.

                                                    {  /*{{{*/
 
   /*Geometry core where we compute indices into tables pre computed in the SealevelRiseAnalysis: */
 
   /*parameters: */
   bool spherical=true;
   IssmDouble *latitude  = NULL;
   IssmDouble *longitude = NULL;
   IssmDouble *radius    = NULL;
   IssmDouble *xx    = NULL;
   IssmDouble *yy    = NULL;
   IssmDouble *zz    = NULL;
   int  horiz;
   bool geometrydone = false;
 
      
   /*retrieve parameters:*/
   femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
   femmodel->parameters->FindParam(&geometrydone,SealevelriseGeometryDoneEnum);
 
   if(geometrydone){
      if(VerboseSolution()) _printf0_("     geometrical offsets have already been computed, skipping \n");
      return; //don't need to run this again.
   }
 
   /*Verbose: */
   if(VerboseSolution()) _printf0_("     computing geometrical offsets into precomputed Green tables \n");
 
   /*first, recover lat,long and radius vectors from vertices: */
   VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); 
   if(horiz) VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); 
 
 
   /*Run sealevelrie geometry routine in elements:*/
   for(int i=0;i<femmodel->elements->Size();i++){
      Element*   element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      element->SealevelriseGeometry(latitude,longitude,radius,xx,yy,zz);
   }
 
   /*Free ressources:*/
   if(horiz){
      xDelete<IssmDouble>(xx);
      xDelete<IssmDouble>(yy);
      xDelete<IssmDouble>(zz);
   }
   xDelete<IssmDouble>(longitude);
   xDelete<IssmDouble>(radius);
 
   /*Record the fact that we ran this module already: */
   femmodel->parameters->SetParam(true,SealevelriseGeometryDoneEnum); 
 
 
}/*}}}*/

sealevelrise_core_eustatic()

Vector<IssmDouble>* sealevelrise_core_eustatic	(	FemModel *	femmodel,
		SealevelMasks *	masks,
		IssmDouble *	poceanarea
	)

Definition at line 404 of file sealevelchange_core.cpp.

                                                                                                               { /*{{{*/
 
   /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/
 
   Vector<IssmDouble> *RSLgi    = NULL;
   IssmDouble          RSLgi_oceanaverage   = 0;
 
   /*parameters: */
   int  gsize;
   IssmDouble oceanarea;
 
   /*outputs:*/
   IssmDouble eustatic;
   
   if(VerboseSolution()) _printf0_("     computing eustatic components on ice\n");
 
   
   /*Figure out size of g-set deflection vector and allocate solution vector: */
   gsize = femmodel->nodes->NumberOfDofs(GsetEnum);
 
   /*Initialize:*/
   RSLgi = new Vector<IssmDouble>(gsize);
 
   /*call the eustatic main module: */
   femmodel->SealevelriseEustatic(RSLgi,&oceanarea,&eustatic, masks); //this computes 
 
   /*we need to average RSLgi over the ocean: RHS term  4 in Eq.4 of Farrel and clarke. Only the elements can do that: */
   RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi,masks, oceanarea);
 
   /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the 
    * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/
   RSLgi->Shift(-eustatic-RSLgi_oceanaverage);
 
   /*save eustatic value for results: */
   femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelRSLEustaticEnum,-eustatic));
 
   /*Assign output pointers and return: */
   *poceanarea=oceanarea;
   return RSLgi;
}/*}}}*/

sealevelrise_core_noneustatic()

Vector<IssmDouble>* sealevelrise_core_noneustatic	(	FemModel *	femmodel,
		SealevelMasks *	masks,
		Vector< IssmDouble > *	RSLg_eustatic,
		IssmDouble	oceanarea
	)

Definition at line 444 of file sealevelchange_core.cpp.

                                                                                                                                                   { /*{{{*/
 
   /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
     non eustatic core of the SLR solution */
 
   Vector<IssmDouble> *RSLg    = NULL;
   Vector<IssmDouble> *RSLg_old    = NULL;
   Vector<IssmDouble> *BPg    = NULL;
 
   Vector<IssmDouble> *RSLgo    = NULL; //ocean convolution of the perturbation to gravity potential.
   Vector<IssmDouble> *RSLgo_rot= NULL; // rotational feedback 
   IssmDouble          RSLgo_oceanaverage = 0;  //average of RSLgo over the ocean.
 
   /*parameters: */
   int count;
   bool save_results;
   int  gsize;
   bool converged=true;
   bool rotation=true;
   bool verboseconvolution=true;
   int max_nonlinear_iterations;
   IssmDouble           eps_rel;
   IssmDouble           eps_abs;
   IssmDouble           eustatic;
   IssmDouble        Ixz, Iyz, Izz; 
   int bp_compute_fingerprints= 0;
   
   if(VerboseSolution()) _printf0_("     converging on ocean components\n");
 
   /*Recover some parameters: */
   femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum);
   femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
   femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
 
   /*computational flag: */
   femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
 
   /*Figure out size of g-set deflection vector and allocate solution vector: */
   gsize = femmodel->nodes->NumberOfDofs(GsetEnum);
 
   /*Initialize:*/
   RSLg = new Vector<IssmDouble>(gsize);
   RSLg->Assemble();
   RSLg_eustatic->Copy(RSLg);  //first initialize RSLg with the eustatic component computed in sealevelrise_core_eustatic.
 
   RSLg_old = new Vector<IssmDouble>(gsize);
   RSLg_old->Assemble();
 
   count=1;
   converged=false;
 
   /*Start loop: */
   for(;;){
 
      //save pointer to old sea level rise
      delete RSLg_old; RSLg_old=RSLg; 
 
      /*Initialize solution vector: */
      RSLg  = new Vector<IssmDouble>(gsize); RSLg->Assemble();
      RSLgo = new Vector<IssmDouble>(gsize); RSLgo->Assemble();
 
      /*call the non eustatic module: */
      femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old,  masks, verboseconvolution);
 
      /*assemble solution vector: */
      RSLgo->Assemble(); 
 
      if(rotation){
 
         /*call rotational feedback  module: */
         RSLgo_rot = new Vector<IssmDouble>(gsize); RSLgo_rot->Assemble();
         femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, masks); 
         RSLgo_rot->Assemble(); 
 
         /*save changes in inertia tensor as results: */
         femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz));
         femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,Iyz));
         femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorZZEnum,Izz));
 
         RSLgo->AXPY(RSLgo_rot,1); 
      }
 
      /*we need to average RSLgo over the ocean: RHS term  5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
      RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo,masks, oceanarea);
 
      /*RSLg is the sum of the eustatic term, and the ocean terms: */
      RSLg_eustatic->Copy(RSLg); RSLg->AXPY(RSLgo,1); 
      RSLg->Shift(-RSLgo_oceanaverage);
 
      /*convergence criterion:*/
      slrconvergence(&converged,RSLg,RSLg_old,eps_rel,eps_abs);
      
   
      /*free ressources: */
      delete RSLgo;
      delete RSLgo_rot;
 
      /*Increase count: */
      count++;
      if(converged==true){
         break;
      }
      if(count>=max_nonlinear_iterations){
         _printf0_("   maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); 
         converged=true;
         break;
      }  
 
      /*some minor verbosing adjustment:*/
      if(count>1)verboseconvolution=false;
 
   }
   if(VerboseConvergence()) _printf0_("\n              total number of iterations: " << count-1 << "\n");
   
   
   /*if we had bottom pressure loading, add dynamic sea level 
    * to RSL:*/
   femmodel->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
   if(bp_compute_fingerprints){
      GetVectorFromInputsx(&BPg,femmodel,DslSeaWaterPressureChangeAtSeaFloor,VertexSIdEnum);
      RSLg->AXPY(BPg,1);
   }
 
   delete RSLg_old;
   delete BPg;
 
   return RSLg;
} /*}}}*/

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 at line 572 of file sealevelchange_core.cpp.

                                                                                                                                                                                                 { /*{{{*/
 
   Vector<IssmDouble> *U_esa  = NULL; 
   Vector<IssmDouble> *U_north_esa   = NULL; 
   Vector<IssmDouble> *U_east_esa    = NULL; 
 
   /*parameters: */
   int  gsize;
   bool spherical=true;
 
   IssmDouble          *latitude   = NULL;
   IssmDouble          *longitude  = NULL;
   IssmDouble          *radius     = NULL;
   IssmDouble          *xx     = NULL;
   IssmDouble          *yy     = NULL;
   IssmDouble          *zz     = NULL;
   int  horiz;
   
   if(VerboseSolution()) _printf0_("     computing vertical and horizontal geodetic signatures\n");
 
   /*retrieve some parameters:*/
   femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
 
   /*find size of vectors:*/
   gsize      = femmodel->nodes->NumberOfDofs(GsetEnum);
 
   /*intialize vectors:*/
   U_esa = new Vector<IssmDouble>(gsize);
   if (horiz){
      U_north_esa = new Vector<IssmDouble>(gsize);
      U_east_esa = new Vector<IssmDouble>(gsize);
   }
 
   /*retrieve geometric information: */
   VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); 
   VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); 
 
   /*call the elastic main modlule:*/ 
   femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg, masks);
 
   /*Assign output pointers:*/
   *pU_esa=U_esa;
   if(horiz){
      *pU_east_esa=U_east_esa;
      *pU_north_esa=U_north_esa;
   }
 
   /*Free ressources: */
   xDelete<IssmDouble>(longitude);
   xDelete<IssmDouble>(latitude);
   xDelete<IssmDouble>(xx);
   xDelete<IssmDouble>(yy);
   xDelete<IssmDouble>(zz);
   xDelete<IssmDouble>(radius);
}

sealevelrise_core_viscous()

void sealevelrise_core_viscous	(	Vector< IssmDouble > **	pU_gia,
		Vector< IssmDouble > **	pN_gia,
		FemModel *	femmodel,
		Vector< IssmDouble > *	RSLg
	)

Definition at line 628 of file sealevelchange_core.cpp.

                                                                                                                                    { /*{{{*/
 
   /*variables:*/
   Vector<IssmDouble> *U_gia  = NULL; 
   Vector<IssmDouble> *N_gia  = NULL; 
 
   /*parameters:*/
   int               frequency;
   IssmDouble          dt;
   
   if(VerboseSolution()) _printf0_("     computing viscous components\n");
 
   /*retrieve some parameters:*/
   femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
 
   /*recover GIA rates:*/
   GetVectorFromInputsx(&U_gia,femmodel,UGiaRateEnum,VertexSIdEnum);
   GetVectorFromInputsx(&N_gia,femmodel,NGiaRateEnum,VertexSIdEnum);
 
   /*we just loaded rates, that's not what's being asked, scale by time:*/ 
   U_gia->Scale(frequency*dt);
   N_gia->Scale(frequency*dt);
 
   /*Assign output pointers:*/
   *pU_gia=U_gia;
   *pN_gia=N_gia;
 
}

sealevelrise_diagnostics()

void sealevelrise_diagnostics	(	FemModel *	femmodel,
		Vector< IssmDouble > *	RSLg
	)

Definition at line 658 of file sealevelchange_core.cpp.

                                                                          { /*{{{*/
 
   /*compute spherical harmonics deg 1 and deg 2 coefficeints:*/
 
}

GetDynamicRate()

void GetDynamicRate	(	Vector< IssmDouble > **	pdynamic_rate_g,
		FemModel *	femmodel
	)

Definition at line 664 of file sealevelchange_core.cpp.

                                                                              { /*{{{*/
 
   int dslmodel=-1;
   IssmDouble time;
 
   /*variables:*/
   Vector<IssmDouble> *dynamic_rate_g  = NULL; 
 
   /*Update steric rates before retrieving them on Vertex SID set:*/
   femmodel->parameters->FindParam(&dslmodel,DslModelEnum);
   femmodel->parameters->FindParam(&time,TimeEnum);
   if(dslmodel==1){
      TransientInput2* transient_input  = femmodel->inputs2->GetTransientInput(DslSeaSurfaceHeightChangeAboveGeoidEnum);
      TriaInput2* tria_input=transient_input->GetTriaInput(time);
      Input2* tria_input_copy=tria_input->copy();
      tria_input_copy->ChangeEnum(DslDynamicRateEnum);
      femmodel->inputs2->AddInput(tria_input_copy);
   }
   else if(dslmodel==2){
   
      int modelid;
      
      /*Recover modelid:*/
      femmodel->parameters->FindParam(&modelid,DslModelidEnum);
      modelid--; //from matlab. 
      
      /*find the DslSeaSurfaceHeightChangeAboveGeoidEnum dataset of transient inputs:*/
      DatasetInput2* dataset_input=femmodel->inputs2->GetDatasetInput2(DslSeaSurfaceHeightChangeAboveGeoidEnum);
      
      /*Go find the modelid'th transient input:*/
      TriaInput2* tria_input=dataset_input->GetTriaInputByOffset(modelid);
      
      /*Plug into DslDynamicRate input: */
      Input2* tria_input_copy=tria_input->copy();
      tria_input_copy->ChangeEnum(DslDynamicRateEnum);
      femmodel->inputs2->AddInput(tria_input_copy);
   }
   else _error_("not implemented yet");
 
   GetVectorFromInputsx(&dynamic_rate_g,femmodel,DslDynamicRateEnum,VertexSIdEnum);
   *pdynamic_rate_g=dynamic_rate_g;
}

GetStericRate()

void GetStericRate	(	Vector< IssmDouble > **	psteric_rate_g,
		FemModel *	femmodel
	)

Definition at line 707 of file sealevelchange_core.cpp.

                                                                            { /*{{{*/
 
   int dslmodel=-1;
   IssmDouble time;
 
   /*variables:*/
   Vector<IssmDouble> *steric_rate_g  = NULL; 
 
   /*Update steric rates before retrieving them on Vertex SID set:*/
   femmodel->parameters->FindParam(&dslmodel,DslModelEnum);
   femmodel->parameters->FindParam(&time,TimeEnum);
   if(dslmodel==1){
      TransientInput2* transient_input  = femmodel->inputs2->GetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum);
      TriaInput2* tria_input=transient_input->GetTriaInput(time);
      Input2* tria_input_copy=tria_input->copy();
      tria_input_copy->ChangeEnum(DslStericRateEnum);
      femmodel->inputs2->AddInput(tria_input_copy);
   }
   else if (dslmodel==2){
      int modelid;
      
      /*Recover modelid:*/
      femmodel->parameters->FindParam(&modelid,DslModelidEnum);
      
      modelid--; //from matlab. 
      
      /*find the DslGlobalAverageThermostericSeaLevelChangeEnum dataset of transient inputs:*/
      DatasetInput2* dataset_input=femmodel->inputs2->GetDatasetInput2(DslGlobalAverageThermostericSeaLevelChangeEnum);
      
      /*Go find the modelid'th transient input:*/
      TriaInput2* tria_input=dataset_input->GetTriaInputByOffset(modelid);
      
      /*Plug into DslStericRate input: */
      Input2* tria_input_copy=tria_input->copy();
      tria_input_copy->ChangeEnum(DslStericRateEnum);
      femmodel->inputs2->AddInput(tria_input_copy);
   }
   else _error_("not implemented yet");
 
   GetVectorFromInputsx(&steric_rate_g,femmodel,DslStericRateEnum,VertexSIdEnum);
   *psteric_rate_g=steric_rate_g;
}

TransferForcing()

void TransferForcing	(	FemModel *	femmodel,
		int	forcingenum
	)

Definition at line 752 of file sealevelchange_core.cpp.

                                                        { /*{{{*/
 
   /*forcing being transferred from models to earth: */
   IssmDouble** forcings=NULL;
   IssmDouble*  forcing=NULL; 
   Vector<IssmDouble>* forcingglobal=NULL; 
   int*         nvs=NULL;
 
   /*transition vectors:*/
   IssmDouble** transitions=NULL;
   int          ntransitions; 
   int*         transitions_m=NULL;
   int*         transitions_n=NULL;
   int          nv;
 
   /*communicators:*/
   ISSM_MPI_Comm tocomm;
   ISSM_MPI_Comm* fromcomms=NULL;
   ISSM_MPI_Status status;
   int         my_rank;
   int         modelid,earthid;
   int         nummodels;
 
   /*Recover some parameters: */
   femmodel->parameters->FindParam(&modelid,ModelIdEnum);
   femmodel->parameters->FindParam(&earthid,EarthIdEnum);
   femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
   my_rank=IssmComm::GetRank();
 
   /*retrieve the inter communicators that will be used to send data from each ice cap to the earth: */
   if(modelid==earthid){
      GenericParam<ISSM_MPI_Comm*>* parcoms = dynamic_cast<GenericParam<ISSM_MPI_Comm*>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
      if(!parcoms)_error_("TransferForcing error message: could not find IcecapToEarthComm communicator");
      fromcomms=parcoms->GetParameterValue();
   }
   else {
      GenericParam<ISSM_MPI_Comm>* parcom = dynamic_cast<GenericParam<ISSM_MPI_Comm>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
      if(!parcom)_error_("TransferForcing error message: could not find IcecapToEarthComm communicator");
      tocomm=parcom->GetParameterValue();
   }
 
   /*For each icecap, retrieve the forcing vector that will be sent to the earth model: */
   if(modelid!=earthid){
      nv=femmodel->vertices->NumberOfVertices();
      GetVectorFromInputsx(&forcing,femmodel,forcingenum,VertexSIdEnum);
   }
 
   /*Send the forcing to the earth model:{{{*/
   if(my_rank==0){
      if(modelid==earthid){
         forcings=xNew<IssmDouble*>(nummodels-1);
         nvs=xNew<int>(nummodels-1);
         for(int i=0;i<earthid;i++){
            ISSM_MPI_Recv(nvs+i, 1, ISSM_MPI_INT, 0,i, fromcomms[i], &status);
            forcings[i]=xNew<IssmDouble>(nvs[i]);
            ISSM_MPI_Recv(forcings[i], nvs[i], ISSM_MPI_DOUBLE, 0,i, fromcomms[i], &status);
         }
 
      }
      else{
         ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, modelid, tocomm);
         ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, tocomm);
      }
   }
   /*}}}*/
 
   /*On the earth model, consolidate all the forcings into one, and update the elements dataset accordingly: {{{*/
   if(modelid==earthid){
 
      /*Out of all the delta thicknesses, build one delta thickness vector made of all the ice cap contributions. 
       *First, build the global delta thickness vector in the earth model: */
      nv=femmodel->vertices->NumberOfVertices();
      GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum);
 
      /*Retrieve transition vectors, used to plug from each ice cap into the global forcing:*/
      femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);
 
      if(ntransitions!=earthid)_error_("TransferForcing error message: number of transition vectors is not equal to the number of icecaps!");
 
      /*Go through all the delta thicknesses coming from each ice cap: */
      if(my_rank==0){
         for(int i=0;i<earthid;i++){
 
            IssmDouble* forcingfromcap= forcings[i]; //careful, this only exists on rank 0 of the earth model!
            IssmDouble* transition=transitions[i];
            int         M=transitions_m[i];
 
            /*build index to plug values: */
            int*        index=xNew<int>(M); for(int i=0;i<M;i++)index[i]=reCast<int>(transition[i])-1; //matlab indexing!
 
            /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular 
             * ice cap: */
            forcingglobal->SetValues(M,index,forcingfromcap,ADD_VAL);
            xDelete<int>(index);
         }
      }
 
      /*Assemble vector:*/
      forcingglobal->Assemble();
 
      /*Plug into elements:*/
      InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum);
   } 
   /*}}}*/
 
   /*Free ressources:{{{*/
   if(forcings){
      for(int i=0;i<nummodels-1;i++){
         IssmDouble* temp=forcings[i]; xDelete<IssmDouble>(temp);
      }
      xDelete<IssmDouble*>(forcings);
   }
   if(forcing)xDelete<IssmDouble>(forcing);
   if(forcingglobal)delete forcingglobal;
   if(transitions){
      for(int i=0;i<earthid;i++){
         IssmDouble* temp=transitions[i];
         xDelete<IssmDouble>(temp);
      }
      xDelete<IssmDouble*>(transitions);
      xDelete<int>(transitions_m);
      xDelete<int>(transitions_n);
   }
   if(nvs)xDelete<int>(nvs);
   /*}}}*/
 
} /*}}}*/

TransferSealevel()

void TransferSealevel	(	FemModel *	femmodel,
		int	forcingenum
	)

Definition at line 879 of file sealevelchange_core.cpp.

                                                         { /*{{{*/
 
   /*forcing being transferred from earth to ice caps: */
   IssmDouble*  forcing=NULL; 
   IssmDouble*  forcingglobal=NULL; 
 
   /*transition vectors:*/
   IssmDouble** transitions=NULL;
   int          ntransitions; 
   int*         transitions_m=NULL;
   int*         transitions_n=NULL;
   int          nv;
 
   /*communicators:*/
   ISSM_MPI_Comm fromcomm;
   ISSM_MPI_Comm* tocomms=NULL;
   ISSM_MPI_Status status;
   int         my_rank;
   int         modelid,earthid;
   int         nummodels;
   int         numcoms;
 
   /*Recover some parameters: */
   femmodel->parameters->FindParam(&modelid,ModelIdEnum);
   femmodel->parameters->FindParam(&earthid,EarthIdEnum);
   femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
   my_rank=IssmComm::GetRank();
 
   /*retrieve the inter communicators that will be used to send data from earth to ice caps:*/
   if(modelid==earthid){
      GenericParam<ISSM_MPI_Comm*>* parcoms = dynamic_cast<GenericParam<ISSM_MPI_Comm*>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
      if(!parcoms)_error_("TransferSealevel error message: could not find IcecapToEarthComm communicator");
      tocomms=parcoms->GetParameterValue();
      //femmodel->parameters->FindParam((int**)(&tocomms),&numcoms,IcecapToEarthCommEnum);
   }
   else{
      GenericParam<ISSM_MPI_Comm>* parcom = dynamic_cast<GenericParam<ISSM_MPI_Comm>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
      if(!parcom)_error_("TransferSealevel error message: could not find IcecapToEarthComm communicator");
      fromcomm=parcom->GetParameterValue();
      //femmodel->parameters->FindParam((int*)(&fromcomm), IcecapToEarthCommEnum);
   }
 
   /*Retrieve sea-level on earth model: */
   if(modelid==earthid){
      nv=femmodel->vertices->NumberOfVertices();
      GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum);
   }
 
   /*Send the forcing to the ice caps:{{{*/
   if(my_rank==0){
 
      if(modelid==earthid){
 
         /*Retrieve transition vectors, used to figure out global forcing contribution to each ice cap's own elements: */
         femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);
 
         if(ntransitions!=earthid)_error_("TransferSeaLevel error message: number of transition vectors is not equal to the number of icecaps!");
 
         for(int i=0;i<earthid;i++){
            nv=transitions_m[i];
            forcing=xNew<IssmDouble>(nv);
            IssmDouble* transition=transitions[i];
            for(int j=0;j<nv;j++){
               forcing[j]=forcingglobal[reCast<int>(transition[j])-1];
            }
            ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, i, tocomms[i]);
            ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, i, tocomms[i]);
         }
      }
      else{
         ISSM_MPI_Recv(&nv, 1, ISSM_MPI_INT, 0, modelid, fromcomm, &status);
         forcing=xNew<IssmDouble>(nv);
         ISSM_MPI_Recv(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, fromcomm, &status);
      }
   }
   /*}}}*/
 
   /*On each ice cap, spread the forcing across cpus, and update the elements dataset accordingly: {{{*/
   if(modelid!=earthid){
 
      ISSM_MPI_Bcast(&nv,1,ISSM_MPI_INT,0,IssmComm::GetComm());
      if(my_rank!=0)forcing=xNew<IssmDouble>(nv);
      ISSM_MPI_Bcast(forcing,nv,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
      /*Plug into elements:*/
      InputUpdateFromVectorx(femmodel,forcing,forcingenum,VertexSIdEnum);
   } 
   /*}}}*/
 
   /*Free ressources:{{{*/
   if(forcingglobal)xDelete<IssmDouble>(forcingglobal);
   if(forcing)xDelete<IssmDouble>(forcing);
   if(transitions){
      for(int i=0;i<ntransitions;i++){
         IssmDouble* temp=transitions[i];
         xDelete<IssmDouble>(temp);
      }
      xDelete<IssmDouble*>(transitions);
      xDelete<int>(transitions_m);
      xDelete<int>(transitions_n);
   }
   /*}}}*/
 
} /*}}}*/

EarthMassTransport()

void EarthMassTransport

(

FemModel *

femmodel

)

Definition at line 983 of file sealevelchange_core.cpp.

                                           { /*{{{*/
 
   IssmDouble time,dt;
   Vector<IssmDouble> *oldthickness    = NULL;
   Vector<IssmDouble> *newthickness    = NULL;
   Vector<IssmDouble> *deltathickness    = NULL;
   Vector<IssmDouble> *cumdeltathickness    = NULL;
   int nv;
 
   if(VerboseSolution()) _printf0_("              computing earth mass transport\n");
 
   /*This module has to be recoded! We do not grab spc thicknesses from an earth model 
    * anymore! Thicknesses come from a mass transport module applied to each basin! 
    * Commeting out for now:*/
   _error_("EarthMassTransport error message: not supported anymore!");
 
   /*This mass transport module for the Earth is because we might have thickness variations as spcs 
    * specified in the md.slr class, outside of what we will get from the icecaps. That's why we get t
    * the thickness variations from SealevelriseSpcthicknessEnum.*/
 
   /*No mass transport module was called, so we are just going to retrieve the geometry thickness 
    * at this time step, at prior time step, and plug the difference as deltathickness: */
   /*femmodel->parameters->FindParam(&time,TimeEnum);
   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
   nv=femmodel->vertices->NumberOfVertices();
 
   GetVectorFromInputsx(&newthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum);
   GetVectorFromInputsx(&oldthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum,time-dt);*/
 
   /*compute deltathickness: */
   /*deltathickness = new Vector<IssmDouble>(nv); 
   newthickness->Copy(deltathickness); deltathickness->AXPY(oldthickness,-1); */
 
   /*plug into elements:*/
// InputUpdateFromVectorx(femmodel,deltathickness,SurfaceloadIceThicknessChangeEnum,VertexSIdEnum);
 
   /*add to cumulated delta thickness: */
   /*GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
   cumdeltathickness->AXPY(deltathickness,1); 
   InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);*/
 
   /*free ressources:*/
   /*delete oldthickness;
   delete newthickness;
   delete deltathickness;
   delete cumdeltathickness;*/
 
} /*}}}*/

slrconvergence()

void slrconvergence	(	bool *	pconverged,
		Vector< IssmDouble > *	RSLg,
		Vector< IssmDouble > *	RSLg_old,
		IssmDouble	eps_rel,
		IssmDouble	eps_abs
	)

Definition at line 1031 of file sealevelchange_core.cpp.

                                                                                                                                  { /*{{{*/
 
   bool converged=true;
   IssmDouble ndS,nS; 
   Vector<IssmDouble> *dRSLg    = NULL;
 
   //compute norm(du) and norm(u) if requested
   dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0);
   ndS=dRSLg->Norm(NORM_TWO); 
 
   if (xIsNan<IssmDouble>(ndS)) _error_("convergence criterion is NaN!");
 
   if(!xIsNan<IssmDouble>(eps_rel)){
      nS=RSLg_old->Norm(NORM_TWO);
      if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN!");
   }
 
   //clean up
   delete dRSLg;
 
   //print
   if(!xIsNan<IssmDouble>(eps_rel)){
      if((ndS/nS)<eps_rel){
         if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " < " << eps_rel*100 << " %\n");
      }
      else{ 
         if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " > " << eps_rel*100 << " %\n");
         converged=false;
      }
   }
   if(!xIsNan<IssmDouble>(eps_abs)){
      if(ndS<eps_abs){
         if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)" << ndS << " < " << eps_abs << " \n");
      }
      else{ 
         if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)" << ndS << " > " << eps_abs << " \n");
         converged=false;
      }
   }
 
   /*assign output*/
   *pconverged=converged;
 
} /*}}}*/