Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/Makefile.am =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/Makefile.am (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/Makefile.am (revision 22155) @@ -476,6 +476,4 @@ if SEALEVELRISE issm_sources += ./cores/sealevelrise_core.cpp\ - ./cores/sealevelrise_core_eustatic.cpp\ - ./cores/sealevelrise_core_noneustatic.cpp\ ./analyses/SealevelriseAnalysis.cpp endif Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp (revision 22155) @@ -20,4 +20,5 @@ void SealevelriseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + int geodetic=0; /*Update elements: */ @@ -34,9 +35,14 @@ iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum); iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); - iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum); - iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum); + //those only if we have requested geodetic computations: + iomodel->FetchData(&geodetic,"md.slr.geodetic"); + if (geodetic){ + iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum); + iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum); + } iomodel->FetchDataToInput(elements,"md.slr.deltathickness",SealevelriseDeltathicknessEnum); iomodel->FetchDataToInput(elements,"md.slr.spcthickness",SealevelriseSpcthicknessEnum); iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(elements,"md.geometry.bed",BedEnum); /*Initialize cumdeltalthickness input: unfortunately, we don't have femmodel, so we need to iterate on @@ -96,4 +102,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic",SealevelriseGeodeticEnum)); iomodel->FetchData(&elastic,"md.slr.elastic"); Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.cpp (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.cpp (revision 22155) @@ -1642,4 +1642,5 @@ name==SealevelEnum || name==SealevelUmotionEnum || + name==SealevelUmotionRateEnum || name==SealevelNmotionEnum || name==SealevelEmotionEnum || @@ -1648,4 +1649,5 @@ name==SealevelriseDeltathicknessEnum || name==SealevelriseCumDeltathicknessEnum || + name==SealevelRateEnum || name==EsaUmotionEnum || name==EsaNmotionEnum || Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h (revision 22155) @@ -51,6 +51,9 @@ void damage_core(FemModel* femmodel); void sealevelrise_core(FemModel* femmodel); +void geodetic_core(FemModel* femmodel); +void steric_core(FemModel* femmodel); Vector* sealevelrise_core_eustatic(FemModel* femmodel); Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* Sg_eustatic); +void sealevelrise_core_elastic(Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* Sg); IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel); @@ -66,4 +69,5 @@ void TransferSealevel(FemModel* femmodel,int forcingenum); void EarthMassTransport(FemModel* femmodel); +void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs); //solution configuration Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp (revision 22155) @@ -10,23 +10,417 @@ #include "../solutionsequences/solutionsequences.h" +/*cores:*/ void sealevelrise_core(FemModel* femmodel){ /*{{{*/ + + /*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 geodetic:*/ + geodetic_core(femmodel); + + /*Run steric core for sure:*/ + steric_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(); +} +/*}}}*/ +void geodetic_core(FemModel* femmodel){ /*{{{*/ + + + /*variables:*/ Vector *Sg = NULL; + Vector *Sg_rate = NULL; + Vector *Sg_eustatic = NULL; + Vector *U_radial = NULL; + Vector *U_radial_rate = NULL; + Vector *U_north = NULL; + Vector *U_east = NULL; + Vector* cumdeltathickness=NULL; + + /*parameters:*/ + bool iscoupler; + int configuration_type; + int modelid,earthid; + bool istransientmasstransport; + int frequency,count; + int horiz; + int geodetic=0; + IssmDouble dt; + + /*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(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + + 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 eustatic core (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */ + Sg_eustatic=sealevelrise_core_eustatic(femmodel); + + /*call non-eustatic core (ocean loading tems - 2nd and 5th terms on the RHS of Farrel and Clark) */ + Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); + + + /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */ + sealevelrise_core_elastic(&U_radial,&U_north,&U_east,femmodel,Sg); + + /*transform these values into rates (as we only run this once each frequency turns:*/ + Sg_rate=Sg->Duplicate(); Sg->Copy(Sg_rate); Sg_rate->Scale(1/(dt*frequency)); + U_radial_rate=U_radial->Duplicate(); U_radial->Copy(U_radial_rate); U_radial_rate->Scale(1/(dt*frequency)); + + + /*get some results into elements:*/ + InputUpdateFromVectorx(femmodel,Sg_rate,SealevelRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,U_radial_rate,SealevelUmotionRateEnum,VertexSIdEnum); + if (horiz){ + InputUpdateFromVectorx(femmodel,U_north,SealevelNmotionEnum,VertexSIdEnum); // north motion + InputUpdateFromVectorx(femmodel,U_east,SealevelEmotionEnum,VertexSIdEnum); // east motion + } + } + + + if(iscoupler){ + /*transfer sea level back to ice caps:*/ + TransferSealevel(femmodel,SealevelRateEnum); + TransferSealevel(femmodel,SealevelUmotionRateEnum); + + //reset cumdeltathickness to 0: + InputUpdateFromConstantx(femmodel,0,SealevelriseCumDeltathicknessEnum); + } + + /*reset counter to 1:*/ + femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter. + + /*Free ressources:*/ + delete Sg_eustatic; + delete Sg; + delete Sg_rate; + delete U_radial; + delete U_radial_rate; + if(horiz){ + delete U_north; + delete U_east; + } +} +/*}}}*/ +void steric_core(FemModel* femmodel){ /*{{{*/ + + /*variables:*/ + Vector *bedrock = NULL; Vector *Sg_absolute = NULL; - Vector *Sg_eustatic = NULL; Vector *steric_rate_g = NULL; + Vector *rsl_g = NULL; + Vector *bedrock_rate_g = 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 absolute sea level, RSL rate, bedrock uplift rate and steric rate, as vectors:*/ + GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum); + GetVectorFromInputsx(&Sg_absolute,femmodel,SealevelEnum,VertexSIdEnum); + GetVectorFromInputsx(&steric_rate_g,femmodel,SealevelriseStericRateEnum,VertexSIdEnum); + if(geodetic){ + GetVectorFromInputsx(&rsl_g,femmodel,SealevelRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&bedrock_rate_g,femmodel,SealevelUmotionRateEnum,VertexSIdEnum); + } + + /*compute: absolute sea level change = initial absolute sea level + relative sea level change rate * dt + vertical motion rate * dt + steric_rate * dt*/ + if(rsl_g) Sg_absolute->AXPY(rsl_g,dt); + if(bedrock_rate_g)Sg_absolute->AXPY(bedrock_rate_g,dt); + Sg_absolute->AXPY(steric_rate_g,dt); + + /*compute new bedrock position: */ + if(bedrock_rate_g)bedrock->AXPY(bedrock_rate_g,dt); + + /*update element inputs:*/ + InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,bedrock,SealevelUmotionEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelEnum,VertexSIdEnum); + + /*Free ressources:*/ + delete bedrock; + delete Sg_absolute; + delete steric_rate_g; + if(rsl_g)delete rsl_g; + if(bedrock_rate_g)delete bedrock_rate_g; + +} +/*}}}*/ +Vector* sealevelrise_core_eustatic(FemModel* femmodel){ /*{{{*/ + + /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/ + + Vector *Sgi = NULL; + IssmDouble Sgi_oceanaverage = 0; + + /*parameters: */ + int configuration_type; + int gsize; + bool spherical=true; + IssmDouble *latitude = NULL; + IssmDouble *longitude = NULL; + IssmDouble *radius = NULL; + int loop; + + /*outputs:*/ + IssmDouble eustatic; + + /*recover parameters:*/ + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum); + + /*first, recover lat,long and radius vectors from vertices: */ + VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); + + /*Figure out size of g-set deflection vector and allocate solution vector: */ + gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); + + /*Initialize:*/ + Sgi = new Vector(gsize); + + /*call the eustatic main module: */ + femmodel->SealevelriseEustatic(Sgi,&eustatic, latitude, longitude, radius,loop); //this computes + + /*we need to average Sgi over the ocean: RHS term 4 in Eq.4 of Farrel and clarke. Only the elements can do that: */ + Sgi_oceanaverage=femmodel->SealevelriseOceanAverage(Sgi); + + /*Sg 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):*/ + Sgi->Shift(-eustatic-Sgi_oceanaverage); + + /*save eustatic value for results: */ + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelEustaticEnum,-eustatic)); + + /*clean up and return:*/ + xDelete(latitude); + xDelete(longitude); + xDelete(radius); + return Sgi; +}/*}}}*/ +Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* Sg_eustatic){ /*{{{*/ + + /*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 *Sg = NULL; + Vector *Sg_old = NULL; + + Vector *Sgo = NULL; //ocean convolution of the perturbation to gravity potential. + Vector *Sgo_rot= NULL; // rotational feedback + IssmDouble Sgo_oceanaverage = 0; //average of Sgo over the ocean. + + /*parameters: */ + int count; + bool save_results; + int gsize; + int configuration_type; + bool spherical=true; + bool converged=true; + bool rotation=true; + bool verboseconvolution=true; + int max_nonlinear_iterations; + IssmDouble eps_rel; + IssmDouble eps_abs; + IssmDouble *latitude = NULL; + IssmDouble *longitude = NULL; + IssmDouble *radius = NULL; + IssmDouble eustatic; + int loop; + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&max_nonlinear_iterations,SealevelriseMaxiterEnum); + femmodel->parameters->FindParam(&eps_rel,SealevelriseReltolEnum); + femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum); + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum); + + /*computational flag: */ + femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum); + + /*first, recover lat,long and radius vectors from vertices: */ + VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); + + /*Figure out size of g-set deflection vector and allocate solution vector: */ + gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); + + /*Initialize:*/ + Sg = new Vector(gsize); + Sg->Assemble(); + Sg_eustatic->Copy(Sg); //first initialize Sg with the eustatic component computed in sealevelrise_core_eustatic. + + Sg_old = new Vector(gsize); + Sg_old->Assemble(); + + count=1; + converged=false; + + /*Start loop: */ + for(;;){ + + //save pointer to old sea level rise + delete Sg_old; Sg_old=Sg; + + /*Initialize solution vector: */ + Sg = new Vector(gsize); Sg->Assemble(); + Sgo = new Vector(gsize); Sgo->Assemble(); + + /*call the non eustatic module: */ + femmodel->SealevelriseNonEustatic(Sgo, Sg_old, latitude, longitude, radius,verboseconvolution,loop); + + /*assemble solution vector: */ + Sgo->Assemble(); + + if(rotation){ + /*call rotational feedback module: */ + Sgo_rot = new Vector(gsize); Sgo_rot->Assemble(); + femmodel->SealevelriseRotationalFeedback(Sgo_rot,Sg_old,latitude,longitude,radius); + Sgo_rot->Assemble(); + + Sgo->AXPY(Sgo_rot,1); + } + + /*we need to average Sgo over the ocean: RHS term 5 in Eq.4 of Farrel and clarke. Only the elements can do that: */ + Sgo_oceanaverage=femmodel->SealevelriseOceanAverage(Sgo); + + /*Sg is the sum of the eustatic term, and the ocean terms: */ + Sg_eustatic->Copy(Sg); Sg->AXPY(Sgo,1); + Sg->Shift(-Sgo_oceanaverage); + + /*convergence criterion:*/ + slrconvergence(&converged,Sg,Sg_old,eps_rel,eps_abs); + + /*free ressources: */ + delete Sgo; + + /*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"); + + xDelete(latitude); + xDelete(longitude); + xDelete(radius); + delete Sg_old; + + return Sg; +} /*}}}*/ +void sealevelrise_core_elastic(Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* Sg){ /*{{{*/ + Vector *U_radial = NULL; Vector *U_north = NULL; Vector *U_east = NULL; - Vector* cumdeltathickness=NULL; - bool save_results,isslr,iscoupler; + + /*parameters: */ int configuration_type; - int solution_type; - int numoutputs = 0; - char **requested_outputs = NULL; - - /*additional parameters: */ int gsize; bool spherical=true; + IssmDouble *latitude = NULL; IssmDouble *longitude = NULL; @@ -35,158 +429,36 @@ IssmDouble *yy = NULL; IssmDouble *zz = NULL; - IssmDouble dt,steric_rate; - int frequency,count; int loop; int horiz; - /*Recover some parameters: */ + /*retrieve some parameters:*/ femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); - femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); - femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); - femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); - femmodel->parameters->FindParam(&frequency,SealevelriseRunFrequencyEnum); - femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum); - - - /*first, recover lat,long and radius vectors from vertices: */ + femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum); + femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum); + + /*find size of vectors:*/ + gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); + + /*intialize vectors:*/ + U_radial = new Vector(gsize); + if (horiz){ + U_north = new Vector(gsize); + U_east = new Vector(gsize); + } + + /*retrieve geometric information: */ VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); - - /*recover x,y,z vectors from vertices: */ VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); - - /*Figure out size of g-set deflection vector and allocate solution vector: */ - gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); - - /*several cases here, depending on value of iscoupler and isslr: - solution_type == SealevelriseSolutionEnum) we are running sea level rise core (no coupler) - ( !iscoupler & !isslr) we are not interested in being here :) - ( !iscoupler & isslr) we are running in uncoupled mode - ( iscoupler & isslr) we are running in coupled mode, and better be earth - ( iscoupler & !isslr) we are running in coupled mode, and better be an ice cap - */ - - if(solution_type==SealevelriseSolutionEnum){ - isslr=1; - iscoupler=0; - count=frequency; - } - - /*early return: */ - if( !iscoupler & !isslr) return; //we are not interested in being here :) - - /*In what follows we assume we are all running slr, either in coupled, or uncoupled mode:*/ - if(VerboseSolution()) _printf0_(" computing sea level rise\n"); - - /*set configuration: */ - if(isslr)femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum); - - /*figure out whether deltathickness was computed on the earth (mass transport module should - * have taken care of it for ice caps: */ - if(iscoupler){ - int modelid,earthid; - bool ismasstransport; - femmodel->parameters->FindParam(&modelid,ModelIdEnum); - femmodel->parameters->FindParam(&earthid,EarthIdEnum); - femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum); - if((modelid==earthid) && (ismasstransport==0)) EarthMassTransport(femmodel); - } - - /*transfer cum deltathickness forcing from ice caps to earth model: */ - if(iscoupler & (count==frequency)) TransferForcing(femmodel,SealevelriseCumDeltathicknessEnum); - - /*call sea-level rise sub cores:*/ - if(isslr & (count==frequency)){ - - if(iscoupler){ - /*we have accumulated thicknesses, dump them in deltathcikness: */ - GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseDeltathicknessEnum,VertexSIdEnum); - delete cumdeltathickness; - } - - femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum); - femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum); - - Sg_eustatic=sealevelrise_core_eustatic(femmodel); //generalized eustatic (Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS. - - Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //ocean loading tems (2nd and 5th terms on the RHS of Farrel and Clark) - - /*compute other geodetic signatures, such as absolute sea level chagne, components of 3-D crustal motion: */ - /*Initialize:*/ - U_radial = new Vector(gsize); - if (horiz){ - U_north = new Vector(gsize); - U_east = new Vector(gsize); - } - Sg_absolute = new Vector(gsize); - - /*call the geodetic main modlule:*/ - femmodel->SealevelriseGeodetic(U_radial,U_north,U_east,Sg,latitude,longitude,radius,xx,yy,zz,loop,horiz); - - /*Now deal with steric ocean expansion by just shifting Sg by a spatial rate pattern : */ - femmodel->parameters->FindParam(&frequency,SealevelriseRunFrequencyEnum); - femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - GetVectorFromInputsx(&steric_rate_g,femmodel,SealevelriseStericRateEnum,VertexSIdEnum); - Sg->AXPY(steric_rate_g,dt*frequency); - - /*compute: absolute sea level change = relative sea level change + vertical motion*/ - Sg->Copy(Sg_absolute); Sg_absolute->AXPY(U_radial,1); - - /*get results into elements:*/ - InputUpdateFromVectorx(femmodel,U_radial,SealevelUmotionEnum,VertexSIdEnum); // radial displacement - if (horiz){ - InputUpdateFromVectorx(femmodel,U_north,SealevelNmotionEnum,VertexSIdEnum); // north motion - InputUpdateFromVectorx(femmodel,U_east,SealevelEmotionEnum,VertexSIdEnum); // east motion - } - InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelAbsoluteEnum,VertexSIdEnum); //absolute sea level - InputUpdateFromVectorx(femmodel,Sg,SealevelEnum,VertexSIdEnum); //relative sea level - - if(save_results){ - if(VerboseSolution()) _printf0_(" saving results\n"); - femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum); - femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); - } - - /*requested dependents: */ - if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx(); - - if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} - - /*Free ressources:*/ - delete U_radial; - if(horiz){ - delete U_north; - delete U_east; - } - delete Sg_absolute; - delete Sg_eustatic; - delete Sg; - delete steric_rate_g; - - } - - /*transfer sea-level back to ice caps, reset cum thicknesses to 0 and reset counter: */ - if(iscoupler){ - if (count==frequency) { - //transfer sea level back to ice caps: - TransferSealevel(femmodel,SealevelEnum); - - //reset cumdeltathickness to 0: - GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); - cumdeltathickness->Set(0); - InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); - delete cumdeltathickness; - - //reset counter to 1: - femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter. - } - else{ - //increment counter: - count++; - femmodel->parameters->SetParam(count,SealevelriseRunCountEnum); - } - } - + + /*call the elastic main modlule:*/ + femmodel->SealevelriseGeodetic(U_radial,U_north,U_east,Sg,latitude,longitude,radius,xx,yy,zz,loop,horiz); + + /*Assign output pointers:*/ + *pU_radial=U_radial; + if(horiz){ + *pU_east=U_east; + *pU_north=U_north; + } + /*Free ressources: */ delete longitude; @@ -196,7 +468,8 @@ delete yy; delete zz; - -} +} /*}}}*/ + +/*support routines:*/ void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/ @@ -475,2 +748,47 @@ } /*}}}*/ +void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ + + bool converged=true; + IssmDouble ndS,nS; + Vector *dSg = NULL; + + //compute norm(du) and norm(u) if requested + dSg=Sg_old->Duplicate(); Sg_old->Copy(dSg); dSg->AYPX(Sg,-1.0); + ndS=dSg->Norm(NORM_TWO); + + if (xIsNan(ndS)) _error_("convergence criterion is NaN!"); + + if(!xIsNan(eps_rel)){ + nS=Sg_old->Norm(NORM_TWO); + if (xIsNan(nS)) _error_("convergence criterion is NaN!"); + } + + + //clean up + delete dSg; + + //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/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core_eustatic.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core_eustatic.cpp (revision 22154) +++ (revision ) @@ -1,61 +1,0 @@ -/*!\file: sealevelrise_core_eustatic.cpp - * \brief: eustatic core of the SLR solution (terms that are constant with respect to sea-level) - */ - -#include "./cores.h" -#include "../toolkits/toolkits.h" -#include "../classes/classes.h" -#include "../shared/shared.h" -#include "../modules/modules.h" -#include "../solutionsequences/solutionsequences.h" - -Vector* sealevelrise_core_eustatic(FemModel* femmodel){ - - Vector *Sgi = NULL; - IssmDouble Sgi_oceanaverage = 0; - - /*parameters: */ - int configuration_type; - int gsize; - bool spherical=true; - IssmDouble *latitude = NULL; - IssmDouble *longitude = NULL; - IssmDouble *radius = NULL; - int loop; - - /*outputs:*/ - IssmDouble eustatic; - - /*recover parameters:*/ - femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); - femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum); - - /*first, recover lat,long and radius vectors from vertices: */ - VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); - - /*Figure out size of g-set deflection vector and allocate solution vector: */ - gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); - - /*Initialize:*/ - Sgi = new Vector(gsize); - - /*call the eustatic main module: */ - femmodel->SealevelriseEustatic(Sgi,&eustatic, latitude, longitude, radius,loop); //this computes - - /*we need to average Sgi over the ocean: RHS term 4 in Eq.4 of Farrel and clarke. Only the elements can do that: */ - Sgi_oceanaverage=femmodel->SealevelriseOceanAverage(Sgi); - - /*Sg 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):*/ - Sgi->Shift(-eustatic-Sgi_oceanaverage); - - /*save eustatic value for results: */ - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelEustaticEnum,-eustatic)); - - /*clean up and return:*/ - xDelete(latitude); - xDelete(longitude); - xDelete(radius); - return Sgi; -} - Index: sm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core_noneustatic.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core_noneustatic.cpp (revision 22154) +++ (revision ) @@ -1,176 +1,0 @@ -/*!\file: sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5. - * \brief: non eustatic core of the SLR solution - */ - -#include "./cores.h" -#include "../toolkits/toolkits.h" -#include "../classes/classes.h" -#include "../shared/shared.h" -#include "../modules/modules.h" -#include "../solutionsequences/solutionsequences.h" - -void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs); - -Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* Sg_eustatic){ /*{{{*/ - - Vector *Sg = NULL; - Vector *Sg_old = NULL; - - Vector *Sgo = NULL; //ocean convolution of the perturbation to gravity potential. - Vector *Sgo_rot= NULL; // rotational feedback - IssmDouble Sgo_oceanaverage = 0; //average of Sgo over the ocean. - - /*parameters: */ - int count; - bool save_results; - int gsize; - int configuration_type; - bool spherical=true; - bool converged=true; - bool rotation=true; - bool verboseconvolution=true; - int max_nonlinear_iterations; - IssmDouble eps_rel; - IssmDouble eps_abs; - IssmDouble *latitude = NULL; - IssmDouble *longitude = NULL; - IssmDouble *radius = NULL; - IssmDouble eustatic; - int loop; - - /*Recover some parameters: */ - femmodel->parameters->FindParam(&max_nonlinear_iterations,SealevelriseMaxiterEnum); - femmodel->parameters->FindParam(&eps_rel,SealevelriseReltolEnum); - femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum); - femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); - femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum); - - /*computational flag: */ - femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum); - - /*first, recover lat,long and radius vectors from vertices: */ - VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); - - /*Figure out size of g-set deflection vector and allocate solution vector: */ - gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); - - /*Initialize:*/ - Sg = new Vector(gsize); - Sg->Assemble(); - Sg_eustatic->Copy(Sg); //first initialize Sg with the eustatic component computed in sealevelrise_core_eustatic. - - Sg_old = new Vector(gsize); - Sg_old->Assemble(); - - count=1; - converged=false; - - /*Start loop: */ - for(;;){ - - //save pointer to old sea level rise - delete Sg_old; Sg_old=Sg; - - /*Initialize solution vector: */ - Sg = new Vector(gsize); Sg->Assemble(); - Sgo = new Vector(gsize); Sgo->Assemble(); - - /*call the non eustatic module: */ - femmodel->SealevelriseNonEustatic(Sgo, Sg_old, latitude, longitude, radius,verboseconvolution,loop); - - /*assemble solution vector: */ - Sgo->Assemble(); - - if(rotation){ - /*call rotational feedback module: */ - Sgo_rot = new Vector(gsize); Sgo_rot->Assemble(); - femmodel->SealevelriseRotationalFeedback(Sgo_rot,Sg_old,latitude,longitude,radius); - Sgo_rot->Assemble(); - - Sgo->AXPY(Sgo_rot,1); - } - - /*we need to average Sgo over the ocean: RHS term 5 in Eq.4 of Farrel and clarke. Only the elements can do that: */ - Sgo_oceanaverage=femmodel->SealevelriseOceanAverage(Sgo); - - /*Sg is the sum of the eustatic term, and the ocean terms: */ - Sg_eustatic->Copy(Sg); Sg->AXPY(Sgo,1); - Sg->Shift(-Sgo_oceanaverage); - - /*convergence criterion:*/ - slrconvergence(&converged,Sg,Sg_old,eps_rel,eps_abs); - - /*free ressources: */ - delete Sgo; - - /*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"); - - xDelete(latitude); - xDelete(longitude); - xDelete(radius); - delete Sg_old; - - return Sg; -} /*}}}*/ - -void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ - - bool converged=true; - IssmDouble ndS,nS; - Vector *dSg = NULL; - - //compute norm(du) and norm(u) if requested - dSg=Sg_old->Duplicate(); Sg_old->Copy(dSg); dSg->AYPX(Sg,-1.0); - ndS=dSg->Norm(NORM_TWO); - - if (xIsNan(ndS)) _error_("convergence criterion is NaN!"); - - if(!xIsNan(eps_rel)){ - nS=Sg_old->Norm(NORM_TWO); - if (xIsNan(nS)) _error_("convergence criterion is NaN!"); - } - - - //clean up - delete dSg; - - //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: /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/transient_core.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/transient_core.cpp (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/transient_core.cpp (revision 22155) @@ -21,5 +21,5 @@ /*parameters: */ IssmDouble finaltime,dt,yts; - bool isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,isslr,iscoupler,ismovingfront,isdamageevolution,ishydrology; + bool isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,isslr,ismovingfront,isdamageevolution,ishydrology; bool save_results,dakota_analysis; bool time_adapt; @@ -55,5 +55,4 @@ femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum); femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); - femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum); femmodel->parameters->FindParam(&ismovingfront,TransientIsmovingfrontEnum); @@ -154,5 +153,5 @@ /*Sea level rise: */ - if(isslr || iscoupler) sealevelrise_core(femmodel); + if(isslr) sealevelrise_core(femmodel); /*unload results*/ Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22155) @@ -66,5 +66,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.inversion.type",InversionTypeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.slr.run_frequency",SealevelriseRunFrequencyEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic_run_frequency",SealevelriseGeodeticRunFrequencyEnum)); parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1)); {/*This is specific to ice...*/ Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumDefinitions.h (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumDefinitions.h (revision 22155) @@ -782,5 +782,7 @@ /*Sea Level Rise{{{*/ SealevelEnum, + SealevelRateEnum, SealevelUmotionEnum, + SealevelUmotionRateEnum, SealevelNmotionEnum, SealevelEmotionEnum, @@ -801,5 +803,6 @@ SealevelriseOceanAreaScalingEnum, SealevelriseStericRateEnum, - SealevelriseRunFrequencyEnum, + SealevelriseGeodeticRunFrequencyEnum, + SealevelriseGeodeticEnum, SealevelriseRunCountEnum, SealevelriseRotationEnum, Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumToStringx.cpp (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumToStringx.cpp (revision 22155) @@ -762,5 +762,7 @@ case LevelsetReinitFrequencyEnum : return "LevelsetReinitFrequency"; case SealevelEnum : return "Sealevel"; + case SealevelRateEnum : return "SealevelRate"; case SealevelUmotionEnum : return "SealevelUmotion"; + case SealevelUmotionRateEnum : return "SealevelUmotionRate"; case SealevelNmotionEnum : return "SealevelNmotion"; case SealevelEmotionEnum : return "SealevelEmotion"; @@ -781,5 +783,6 @@ case SealevelriseOceanAreaScalingEnum : return "SealevelriseOceanAreaScaling"; case SealevelriseStericRateEnum : return "SealevelriseStericRate"; - case SealevelriseRunFrequencyEnum : return "SealevelriseRunFrequency"; + case SealevelriseGeodeticRunFrequencyEnum : return "SealevelriseGeodeticRunFrequency"; + case SealevelriseGeodeticEnum : return "SealevelriseGeodetic"; case SealevelriseRunCountEnum : return "SealevelriseRunCount"; case SealevelriseRotationEnum : return "SealevelriseRotation"; Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/StringToEnumx.cpp (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/StringToEnumx.cpp (revision 22155) @@ -780,5 +780,7 @@ else if (strcmp(name,"LevelsetReinitFrequency")==0) return LevelsetReinitFrequencyEnum; else if (strcmp(name,"Sealevel")==0) return SealevelEnum; + else if (strcmp(name,"SealevelRate")==0) return SealevelRateEnum; else if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum; + else if (strcmp(name,"SealevelUmotionRate")==0) return SealevelUmotionRateEnum; else if (strcmp(name,"SealevelNmotion")==0) return SealevelNmotionEnum; else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum; @@ -799,5 +801,6 @@ else if (strcmp(name,"SealevelriseOceanAreaScaling")==0) return SealevelriseOceanAreaScalingEnum; else if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum; - else if (strcmp(name,"SealevelriseRunFrequency")==0) return SealevelriseRunFrequencyEnum; + else if (strcmp(name,"SealevelriseGeodeticRunFrequency")==0) return SealevelriseGeodeticRunFrequencyEnum; + else if (strcmp(name,"SealevelriseGeodetic")==0) return SealevelriseGeodeticEnum; else if (strcmp(name,"SealevelriseRunCount")==0) return SealevelriseRunCountEnum; else if (strcmp(name,"SealevelriseRotation")==0) return SealevelriseRotationEnum; @@ -872,11 +875,11 @@ else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum; else if (strcmp(name,"Masscon")==0) return MassconEnum; - else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; - else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; - else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"VectorParam")==0) return VectorParamEnum; + if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; + else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; + else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum; + else if (strcmp(name,"VectorParam")==0) return VectorParamEnum; else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; else if (strcmp(name,"Segment")==0) return SegmentEnum; @@ -995,11 +998,11 @@ else if (strcmp(name,"P1")==0) return P1Enum; else if (strcmp(name,"P1DG")==0) return P1DGEnum; - else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum; - else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum; - else if (strcmp(name,"P2")==0) return P2Enum; else stage=9; } if(stage==9){ - if (strcmp(name,"P2bubble")==0) return P2bubbleEnum; + if (strcmp(name,"P1bubble")==0) return P1bubbleEnum; + else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum; + else if (strcmp(name,"P2")==0) return P2Enum; + else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum; else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum; else if (strcmp(name,"P2xP1")==0) return P2xP1Enum; Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/pairoptions.m =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/pairoptions.m (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/pairoptions.m (revision 22155) @@ -127,5 +127,5 @@ for i=1:numoptions, if ~self.list{i,3}, - disp(['WARNING: option ' self.list{i,1} ' was not used']) + %disp(['WARNING: option ' self.list{i,1} ' was not used']) end end Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/sealevelmodel.m =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/sealevelmodel.m (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/sealevelmodel.m (revision 22155) @@ -61,5 +61,5 @@ %check that run_frequency is the same everywhere: for i=1:length(slm.icecaps), - if slm.icecaps{i}.slr.run_frequency~=slm.earth.slr.run_frequency, + if slm.icecaps{i}.slr.geodetic_run_frequency~=slm.earth.slr.geodetic_run_frequency, error(sprintf('sealevelmodel checkconsistenty error: icecap model %s should have the same run frequency as earth!',slm.icecaps{i}.miscellaneous.name)); end Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/slr.m =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/slr.m (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/slr.m (revision 22155) @@ -26,5 +26,6 @@ ocean_area_scaling = 0; steric_rate = 0; %rate of ocean expansion from steric effects. - run_frequency = 1; %how many time steps we skip before we run SLR solver during transient + geodetic_run_frequency = 1; %how many time steps we skip before we run the geodetic part of the solver during transient + geodetic = 0; %compute geodetic SLR? (in addition to steric?) degacc = 0; loop_increment = 0; @@ -53,4 +54,5 @@ %computational flags: + self.geodetic=0; self.rigid=1; self.elastic=1; @@ -80,5 +82,5 @@ %how many time steps we skip before we run SLR solver during transient - self.run_frequency=1; + self.geodetic_run_frequency=1; %output default: @@ -110,5 +112,5 @@ md = checkfield(md,'fieldname','slr.abstol','size',[1 1]); md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1); - md = checkfield(md,'fieldname','slr.run_frequency','size',[1 1],'>=',1); + md = checkfield(md,'fieldname','slr.geodetic_run_frequency','size',[1 1],'>=',1); md = checkfield(md,'fieldname','slr.steric_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10); @@ -129,4 +131,19 @@ error('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!'); end + + %check that if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not, + %a coupler to a planet model is provided. + if self.geodetic, + if md.transient.iscoupler, + %we are good; + else + if strcmpi(class(md.mesh),'mesh3dsurface'), + %we are good + else + error('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!'); + end + end + end + end % }}} @@ -154,10 +171,8 @@ fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); fielddisplay(self,'steric_rate','rate of steric ocean expansion (in mm/yr)'); - fielddisplay(self,'run_frequency','how many time steps we skip before we run SLR solver during transient (default: 1)'); + fielddisplay(self,'geodetic_run_frequency','how many time steps we skip before we run SLR solver during transient (default: 1)'); fielddisplay(self,'rigid','rigid earth graviational potential perturbation'); fielddisplay(self,'elastic','elastic earth graviational potential perturbation'); fielddisplay(self,'rotation','earth rotational potential perturbation'); - fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); - fielddisplay(self,'steric_rate','rate of steric ocean expansion (in mm/yr)'); fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions'); fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'); @@ -186,5 +201,5 @@ WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean'); - WriteData(fid,prefix,'object',self,'fieldname','run_frequency','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','geodetic_run_frequency','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','steric_rate','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double'); @@ -192,4 +207,5 @@ WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer'); %process requested outputs @@ -224,5 +240,5 @@ writejsdouble(fid,[modelname '.slr.rotation'],self.rotation); writejsdouble(fid,[modelname '.slr.ocean_area_scaling'],self.ocean_area_scaling); - writejsdouble(fid,[modelname '.slr.run_frequency'],self.run_frequency); + writejsdouble(fid,[modelname '.slr.geodetic_run_frequency'],self.geodetic_run_frequency); writejsdouble(fid,[modelname '.slr.elastic'],self.elastic); writejs1Darray(fid,[modelname '.slr.steric_rate'],self.steric_rate); Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/solve/solveslm.m =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/m/solve/solveslm.m (revision 22154) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/solve/solveslm.m (revision 22155) @@ -34,5 +34,5 @@ options=pairoptions(varargin{:},'solutionstring',solutionstring); -%figure out if the sum of cluster processors requested sums up correctly: +%make sure we request sum of cluster processors totalnp=0; for i=1:length(slm.icecaps), totalnp=totalnp+slm.icecaps{i}.cluster.np; end