Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 24984) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 24985) @@ -545,5 +545,5 @@ if SEALEVELRISE issm_sources += \ - ./cores/sealevelrise_core.cpp \ + ./cores/sealevelchange_core.cpp \ ./analyses/SealevelriseAnalysis.cpp endif Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 24985) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 24985) @@ -0,0 +1,1060 @@ +/*!\file: sealevelchange_core.cpp + * \brief: core of the sea-level change solution + */ + +#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" + + +/*main cores:*/ +void sealevelchange_core(FemModel* femmodel){ /*{{{*/ + + /*Start profiler*/ + femmodel->profiler->Start(SLRCORE); + + /*Parameters, variables:*/ + bool save_results; + bool isslr=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; + + /*Should we be here?:*/ + if(!isslr)return; + + /*Verbose: */ + if(VerboseSolution()) _printf0_(" computing sea level rise\n"); + + /*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(requested_outputs[i]);} xDelete(requested_outputs);} + } + + /*requested dependents: */ + if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx(); + + /*End profiler*/ + femmodel->profiler->Stop(SLRCORE); +} +/*}}}*/ +void grd_core(FemModel* femmodel){ /*{{{*/ + + /*Gravity rotation deformation core GRD: */ + + /*variables:*/ + Vector *RSLg = NULL; + Vector *BPg = NULL; + Vector *RSLg_rate = NULL; + Vector *RSLg_eustatic = NULL; + Vector *U_esa = NULL; + Vector *U_esa_rate = NULL; + Vector *N_esa = NULL; + Vector *N_esa_rate = NULL; + Vector *U_north_esa = NULL; + Vector *U_east_esa = NULL; + Vector *N_gia= NULL; + Vector *U_gia= NULL; + Vector *N_gia_rate= NULL; + Vector *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,SealevelriseGeodeticEnum); 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,SealevelriseGeodeticRunFrequencyEnum); + femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum); + femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum); + 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 (countparameters->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,SealevelriseDeltathicknessEnum); + } + + /*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,SealevelUGiaRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,N_gia_rate,SealevelNGiaRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,RSLg_rate,SealevelRSLRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,U_esa,SealevelUEsaEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,N_esa,SealevelNEsaEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,U_gia,SealevelUGiaEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,N_gia,SealevelNGiaEnum,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,SealevelNGiaRateEnum); + TransferSealevel(femmodel,SealevelUEsaRateEnum); + TransferSealevel(femmodel,SealevelUGiaRateEnum); + + //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; + /*}}}*/ + +} +/*}}}*/ +void dynstr_core(FemModel* femmodel){ /*{{{*/ + + /*variables:*/ + Vector *bedrock = NULL; + Vector *SL = NULL; + Vector *steric_rate_g = NULL; + Vector *dynamic_rate_g = NULL; + Vector *hydro_rate_g = NULL; + Vector *U_esa_rate= NULL; + Vector *N_esa_rate= NULL; + Vector *U_gia_rate= NULL; + Vector *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,SealevelriseGeodeticEnum); + + /*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); + GetVectorFromInputsx(&hydro_rate_g,femmodel,SealevelriseHydroRateEnum,VertexSIdEnum); + if(geodetic){ + GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&U_gia_rate,femmodel,SealevelUGiaRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&N_esa_rate,femmodel,SealevelNEsaRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&N_gia_rate,femmodel,SealevelNGiaRateEnum,VertexSIdEnum); + } + + /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate) * dt + steric_rate + dynamic_rate + hydro_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); + SL->AXPY(hydro_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; + delete hydro_rate_g; + if(geodetic){ + delete U_esa_rate; + delete U_gia_rate; + delete N_esa_rate; + delete N_gia_rate; + } +} +/*}}}*/ + +SealevelMasks* sealevelrise_core_masks(FemModel* femmodel) { /*{{{*/ + + 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;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->SetSealevelMasks(masks); + } + + return masks; +}/*}}}*/ +void sealevelrise_core_geometry(FemModel* femmodel) { /*{{{*/ + + /*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,SealevelriseHorizEnum); + 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;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->SealevelriseGeometry(latitude,longitude,radius,xx,yy,zz); + } + + /*Free ressources:*/ + if(horiz){ + xDelete(xx); + xDelete(yy); + xDelete(zz); + } + xDelete(longitude); + xDelete(radius); + + /*Record the fact that we ran this module already: */ + femmodel->parameters->SetParam(true,SealevelriseGeometryDoneEnum); + + +}/*}}}*/ +Vector* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* poceanarea){ /*{{{*/ + + /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/ + + Vector *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(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(femmodel->results->Size()+1,SealevelRSLEustaticEnum,-eustatic)); + + /*Assign output pointers and return: */ + *poceanarea=oceanarea; + return RSLgi; +}/*}}}*/ +Vector* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector* RSLg_eustatic,IssmDouble oceanarea){ /*{{{*/ + + /*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 *RSLg = NULL; + Vector *RSLg_old = NULL; + Vector *BPg = NULL; + + Vector *RSLgo = NULL; //ocean convolution of the perturbation to gravity potential. + Vector *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,SealevelriseMaxiterEnum); + femmodel->parameters->FindParam(&eps_rel,SealevelriseReltolEnum); + femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum); + + /*computational flag: */ + femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum); + + /*Figure out size of g-set deflection vector and allocate solution vector: */ + gsize = femmodel->nodes->NumberOfDofs(GsetEnum); + + /*Initialize:*/ + RSLg = new Vector(gsize); + RSLg->Assemble(); + RSLg_eustatic->Copy(RSLg); //first initialize RSLg with the eustatic component computed in sealevelrise_core_eustatic. + + RSLg_old = new Vector(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(gsize); RSLg->Assemble(); + RSLgo = new Vector(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(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(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,Iyz)); + femmodel->results->AddResult(new GenericExternalResult(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; +} /*}}}*/ +void sealevelrise_core_elastic(Vector** pU_esa, Vector** pU_north_esa,Vector** pU_east_esa,FemModel* femmodel,Vector* RSLg, SealevelMasks* masks){ /*{{{*/ + + Vector *U_esa = NULL; + Vector *U_north_esa = NULL; + Vector *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,SealevelriseHorizEnum); + + /*find size of vectors:*/ + gsize = femmodel->nodes->NumberOfDofs(GsetEnum); + + /*intialize vectors:*/ + U_esa = new Vector(gsize); + if (horiz){ + U_north_esa = new Vector(gsize); + U_east_esa = new Vector(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(longitude); + xDelete(latitude); + xDelete(xx); + xDelete(yy); + xDelete(zz); + xDelete(radius); +} +/*}}}*/ +void sealevelrise_core_viscous(Vector** pU_gia, Vector** pN_gia,FemModel* femmodel,Vector* RSLg){ /*{{{*/ + + /*variables:*/ + Vector *U_gia = NULL; + Vector *N_gia = NULL; + + /*parameters:*/ + int frequency; + IssmDouble dt; + + if(VerboseSolution()) _printf0_(" computing viscous components\n"); + + /*retrieve some parameters:*/ + femmodel->parameters->FindParam(&frequency,SealevelriseGeodeticRunFrequencyEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + + /*recover GIA rates:*/ + GetVectorFromInputsx(&U_gia,femmodel,SealevelUGiaRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&N_gia,femmodel,SealevelNGiaRateEnum,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; + +} +/*}}}*/ +void sealevelrise_diagnostics(FemModel* femmodel,Vector* RSLg){ /*{{{*/ + + /*compute spherical harmonics deg 1 and deg 2 coefficeints:*/ + +} +/*}}}*/ +void GetDynamicRate(Vector ** pdynamic_rate_g, FemModel* femmodel){ /*{{{*/ + + int dslmodel=-1; + IssmDouble time; + + /*variables:*/ + Vector *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; +} +/*}}}*/ +void GetStericRate(Vector ** psteric_rate_g, FemModel* femmodel){ /*{{{*/ + + int dslmodel=-1; + IssmDouble time; + + /*variables:*/ + Vector *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; +} +/*}}}*/ + +/*support routines:*/ +void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/ + + /*forcing being transferred from models to earth: */ + IssmDouble** forcings=NULL; + IssmDouble* forcing=NULL; + Vector* 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* parcoms = dynamic_cast*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum)); + if(!parcoms)_error_("TransferForcing error message: could not find IcecapToEarthComm communicator"); + fromcomms=parcoms->GetParameterValue(); + } + else { + GenericParam* parcom = dynamic_cast*>(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(nummodels-1); + nvs=xNew(nummodels-1); + for(int i=0;i(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(M); for(int i=0;i(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(index); + } + } + + /*Assemble vector:*/ + forcingglobal->Assemble(); + + /*Plug into elements:*/ + InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum); + } + /*}}}*/ + + /*Free ressources:{{{*/ + if(forcings){ + for(int i=0;i(temp); + } + xDelete(forcings); + } + if(forcing)xDelete(forcing); + if(forcingglobal)delete forcingglobal; + if(transitions){ + for(int i=0;i(temp); + } + xDelete(transitions); + xDelete(transitions_m); + xDelete(transitions_n); + } + if(nvs)xDelete(nvs); + /*}}}*/ + +} /*}}}*/ +void TransferSealevel(FemModel* femmodel,int forcingenum){ /*{{{*/ + + /*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* parcoms = dynamic_cast*>(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* parcom = dynamic_cast*>(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(nv); + IssmDouble* transition=transitions[i]; + for(int j=0;j(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(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(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(forcingglobal); + if(forcing)xDelete(forcing); + if(transitions){ + for(int i=0;i(temp); + } + xDelete(transitions); + xDelete(transitions_m); + xDelete(transitions_n); + } + /*}}}*/ + +} /*}}}*/ +void EarthMassTransport(FemModel* femmodel){ /*{{{*/ + + IssmDouble time,dt; + Vector *oldthickness = NULL; + Vector *newthickness = NULL; + Vector *deltathickness = NULL; + Vector *cumdeltathickness = NULL; + int nv; + + if(VerboseSolution()) _printf0_(" computing earth mass transport\n"); + + /*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(nv); + newthickness->Copy(deltathickness); deltathickness->AXPY(oldthickness,-1); + + /*plug into elements:*/ + InputUpdateFromVectorx(femmodel,deltathickness,SealevelriseDeltathicknessEnum,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; + +} /*}}}*/ +void slrconvergence(bool* pconverged, Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ + + bool converged=true; + IssmDouble ndS,nS; + Vector *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(ndS)) _error_("convergence criterion is NaN!"); + + if(!xIsNan(eps_rel)){ + nS=RSLg_old->Norm(NORM_TWO); + if (xIsNan(nS)) _error_("convergence criterion is NaN!"); + } + + //clean up + delete dRSLg; + + //print + if(!xIsNan(eps_rel)){ + if((ndS/nS) " << eps_rel*100 << " %\n"); + converged=false; + } + } + if(!xIsNan(eps_abs)){ + if(ndS " << eps_abs << " \n"); + converged=false; + } + } + + /*assign output*/ + *pconverged=converged; + +} /*}}}*/ Index: sm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24984) +++ (revision ) @@ -1,1060 +1,0 @@ -/*!\file: sealevelchange_core.cpp - * \brief: core of the sea-level change solution - */ - -#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" - - -/*main cores:*/ -void sealevelchange_core(FemModel* femmodel){ /*{{{*/ - - /*Start profiler*/ - femmodel->profiler->Start(SLRCORE); - - /*Parameters, variables:*/ - bool save_results; - bool isslr=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; - - /*Should we be here?:*/ - if(!isslr)return; - - /*Verbose: */ - if(VerboseSolution()) _printf0_(" computing sea level rise\n"); - - /*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(requested_outputs[i]);} xDelete(requested_outputs);} - } - - /*requested dependents: */ - if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx(); - - /*End profiler*/ - femmodel->profiler->Stop(SLRCORE); -} -/*}}}*/ -void grd_core(FemModel* femmodel){ /*{{{*/ - - /*Gravity rotation deformation core GRD: */ - - /*variables:*/ - Vector *RSLg = NULL; - Vector *BPg = NULL; - Vector *RSLg_rate = NULL; - Vector *RSLg_eustatic = NULL; - Vector *U_esa = NULL; - Vector *U_esa_rate = NULL; - Vector *N_esa = NULL; - Vector *N_esa_rate = NULL; - Vector *U_north_esa = NULL; - Vector *U_east_esa = NULL; - Vector *N_gia= NULL; - Vector *U_gia= NULL; - Vector *N_gia_rate= NULL; - Vector *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,SealevelriseGeodeticEnum); 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,SealevelriseGeodeticRunFrequencyEnum); - femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum); - femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum); - 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 (countparameters->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,SealevelriseDeltathicknessEnum); - } - - /*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,SealevelUGiaRateEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,N_gia_rate,SealevelNGiaRateEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,RSLg_rate,SealevelRSLRateEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,U_esa,SealevelUEsaEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,N_esa,SealevelNEsaEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,U_gia,SealevelUGiaEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,N_gia,SealevelNGiaEnum,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,SealevelNGiaRateEnum); - TransferSealevel(femmodel,SealevelUEsaRateEnum); - TransferSealevel(femmodel,SealevelUGiaRateEnum); - - //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; - /*}}}*/ - -} -/*}}}*/ -void dynstr_core(FemModel* femmodel){ /*{{{*/ - - /*variables:*/ - Vector *bedrock = NULL; - Vector *SL = NULL; - Vector *steric_rate_g = NULL; - Vector *dynamic_rate_g = NULL; - Vector *hydro_rate_g = NULL; - Vector *U_esa_rate= NULL; - Vector *N_esa_rate= NULL; - Vector *U_gia_rate= NULL; - Vector *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,SealevelriseGeodeticEnum); - - /*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); - GetVectorFromInputsx(&hydro_rate_g,femmodel,SealevelriseHydroRateEnum,VertexSIdEnum); - if(geodetic){ - GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum); - GetVectorFromInputsx(&U_gia_rate,femmodel,SealevelUGiaRateEnum,VertexSIdEnum); - GetVectorFromInputsx(&N_esa_rate,femmodel,SealevelNEsaRateEnum,VertexSIdEnum); - GetVectorFromInputsx(&N_gia_rate,femmodel,SealevelNGiaRateEnum,VertexSIdEnum); - } - - /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate) * dt + steric_rate + dynamic_rate + hydro_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); - SL->AXPY(hydro_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; - delete hydro_rate_g; - if(geodetic){ - delete U_esa_rate; - delete U_gia_rate; - delete N_esa_rate; - delete N_gia_rate; - } -} -/*}}}*/ - -SealevelMasks* sealevelrise_core_masks(FemModel* femmodel) { /*{{{*/ - - 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;ielements->Size();i++){ - Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - element->SetSealevelMasks(masks); - } - - return masks; -}/*}}}*/ -void sealevelrise_core_geometry(FemModel* femmodel) { /*{{{*/ - - /*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,SealevelriseHorizEnum); - 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;ielements->Size();i++){ - Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - element->SealevelriseGeometry(latitude,longitude,radius,xx,yy,zz); - } - - /*Free ressources:*/ - if(horiz){ - xDelete(xx); - xDelete(yy); - xDelete(zz); - } - xDelete(longitude); - xDelete(radius); - - /*Record the fact that we ran this module already: */ - femmodel->parameters->SetParam(true,SealevelriseGeometryDoneEnum); - - -}/*}}}*/ -Vector* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* poceanarea){ /*{{{*/ - - /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/ - - Vector *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(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(femmodel->results->Size()+1,SealevelRSLEustaticEnum,-eustatic)); - - /*Assign output pointers and return: */ - *poceanarea=oceanarea; - return RSLgi; -}/*}}}*/ -Vector* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector* RSLg_eustatic,IssmDouble oceanarea){ /*{{{*/ - - /*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 *RSLg = NULL; - Vector *RSLg_old = NULL; - Vector *BPg = NULL; - - Vector *RSLgo = NULL; //ocean convolution of the perturbation to gravity potential. - Vector *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,SealevelriseMaxiterEnum); - femmodel->parameters->FindParam(&eps_rel,SealevelriseReltolEnum); - femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum); - - /*computational flag: */ - femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum); - - /*Figure out size of g-set deflection vector and allocate solution vector: */ - gsize = femmodel->nodes->NumberOfDofs(GsetEnum); - - /*Initialize:*/ - RSLg = new Vector(gsize); - RSLg->Assemble(); - RSLg_eustatic->Copy(RSLg); //first initialize RSLg with the eustatic component computed in sealevelrise_core_eustatic. - - RSLg_old = new Vector(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(gsize); RSLg->Assemble(); - RSLgo = new Vector(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(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(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz)); - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,Iyz)); - femmodel->results->AddResult(new GenericExternalResult(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; -} /*}}}*/ -void sealevelrise_core_elastic(Vector** pU_esa, Vector** pU_north_esa,Vector** pU_east_esa,FemModel* femmodel,Vector* RSLg, SealevelMasks* masks){ /*{{{*/ - - Vector *U_esa = NULL; - Vector *U_north_esa = NULL; - Vector *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,SealevelriseHorizEnum); - - /*find size of vectors:*/ - gsize = femmodel->nodes->NumberOfDofs(GsetEnum); - - /*intialize vectors:*/ - U_esa = new Vector(gsize); - if (horiz){ - U_north_esa = new Vector(gsize); - U_east_esa = new Vector(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(longitude); - xDelete(latitude); - xDelete(xx); - xDelete(yy); - xDelete(zz); - xDelete(radius); -} -/*}}}*/ -void sealevelrise_core_viscous(Vector** pU_gia, Vector** pN_gia,FemModel* femmodel,Vector* RSLg){ /*{{{*/ - - /*variables:*/ - Vector *U_gia = NULL; - Vector *N_gia = NULL; - - /*parameters:*/ - int frequency; - IssmDouble dt; - - if(VerboseSolution()) _printf0_(" computing viscous components\n"); - - /*retrieve some parameters:*/ - femmodel->parameters->FindParam(&frequency,SealevelriseGeodeticRunFrequencyEnum); - femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - - /*recover GIA rates:*/ - GetVectorFromInputsx(&U_gia,femmodel,SealevelUGiaRateEnum,VertexSIdEnum); - GetVectorFromInputsx(&N_gia,femmodel,SealevelNGiaRateEnum,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; - -} -/*}}}*/ -void sealevelrise_diagnostics(FemModel* femmodel,Vector* RSLg){ /*{{{*/ - - /*compute spherical harmonics deg 1 and deg 2 coefficeints:*/ - -} -/*}}}*/ -void GetDynamicRate(Vector ** pdynamic_rate_g, FemModel* femmodel){ /*{{{*/ - - int dslmodel=-1; - IssmDouble time; - - /*variables:*/ - Vector *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; -} -/*}}}*/ -void GetStericRate(Vector ** psteric_rate_g, FemModel* femmodel){ /*{{{*/ - - int dslmodel=-1; - IssmDouble time; - - /*variables:*/ - Vector *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; -} -/*}}}*/ - -/*support routines:*/ -void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/ - - /*forcing being transferred from models to earth: */ - IssmDouble** forcings=NULL; - IssmDouble* forcing=NULL; - Vector* 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* parcoms = dynamic_cast*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum)); - if(!parcoms)_error_("TransferForcing error message: could not find IcecapToEarthComm communicator"); - fromcomms=parcoms->GetParameterValue(); - } - else { - GenericParam* parcom = dynamic_cast*>(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(nummodels-1); - nvs=xNew(nummodels-1); - for(int i=0;i(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(M); for(int i=0;i(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(index); - } - } - - /*Assemble vector:*/ - forcingglobal->Assemble(); - - /*Plug into elements:*/ - InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum); - } - /*}}}*/ - - /*Free ressources:{{{*/ - if(forcings){ - for(int i=0;i(temp); - } - xDelete(forcings); - } - if(forcing)xDelete(forcing); - if(forcingglobal)delete forcingglobal; - if(transitions){ - for(int i=0;i(temp); - } - xDelete(transitions); - xDelete(transitions_m); - xDelete(transitions_n); - } - if(nvs)xDelete(nvs); - /*}}}*/ - -} /*}}}*/ -void TransferSealevel(FemModel* femmodel,int forcingenum){ /*{{{*/ - - /*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* parcoms = dynamic_cast*>(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* parcom = dynamic_cast*>(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(nv); - IssmDouble* transition=transitions[i]; - for(int j=0;j(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(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(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(forcingglobal); - if(forcing)xDelete(forcing); - if(transitions){ - for(int i=0;i(temp); - } - xDelete(transitions); - xDelete(transitions_m); - xDelete(transitions_n); - } - /*}}}*/ - -} /*}}}*/ -void EarthMassTransport(FemModel* femmodel){ /*{{{*/ - - IssmDouble time,dt; - Vector *oldthickness = NULL; - Vector *newthickness = NULL; - Vector *deltathickness = NULL; - Vector *cumdeltathickness = NULL; - int nv; - - if(VerboseSolution()) _printf0_(" computing earth mass transport\n"); - - /*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(nv); - newthickness->Copy(deltathickness); deltathickness->AXPY(oldthickness,-1); - - /*plug into elements:*/ - InputUpdateFromVectorx(femmodel,deltathickness,SealevelriseDeltathicknessEnum,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; - -} /*}}}*/ -void slrconvergence(bool* pconverged, Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ - - bool converged=true; - IssmDouble ndS,nS; - Vector *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(ndS)) _error_("convergence criterion is NaN!"); - - if(!xIsNan(eps_rel)){ - nS=RSLg_old->Norm(NORM_TWO); - if (xIsNan(nS)) _error_("convergence criterion is NaN!"); - } - - //clean up - delete dRSLg; - - //print - if(!xIsNan(eps_rel)){ - if((ndS/nS) " << eps_rel*100 << " %\n"); - converged=false; - } - } - if(!xIsNan(eps_abs)){ - if(ndS " << eps_abs << " \n"); - converged=false; - } - } - - /*assign output*/ - *pconverged=converged; - -} /*}}}*/