Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 25946) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 25947) @@ -552,5 +552,5 @@ issm_sources += \ ./cores/sealevelchange_core.cpp \ - ./analyses/SealevelriseAnalysis.cpp + ./analyses/SealevelchangeAnalysis.cpp endif #}}} Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 25946) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 25947) @@ -105,5 +105,5 @@ #endif #ifdef _HAVE_SEALEVELRISE_ - case SealevelriseAnalysisEnum : return new SealevelriseAnalysis(); + case SealevelchangeAnalysisEnum : return new SealevelchangeAnalysis(); #endif #ifdef _HAVE_SMB_ Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 25946) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 25947) @@ -156,5 +156,5 @@ if(isgroundingline) iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum); /*Initialize cumdeltalthickness input*/ - InputUpdateFromConstantx(inputs,elements,0.,SealevelriseCumDeltathicknessEnum); + InputUpdateFromConstantx(inputs,elements,0.,SealevelchangeCumDeltathicknessEnum); /*Initialize ThicknessResidual input*/ InputUpdateFromConstantx(inputs,elements,0.,ThicknessResidualEnum); @@ -856,5 +856,5 @@ basalelement->GetInputListOnVertices(&phi[0],MaskOceanLevelsetEnum); basalelement->GetInputListOnVertices(&sealevel[0],SealevelEnum); - basalelement->GetInputListOnVertices(&cumdeltathickness[0],SealevelriseCumDeltathicknessOldEnum); + basalelement->GetInputListOnVertices(&cumdeltathickness[0],SealevelchangeCumDeltathicknessOldEnum); /*Do we do grounding line migration?*/ @@ -863,5 +863,5 @@ if(isgroundingline) basalelement->GetInputListOnVertices(&bed[0],BedEnum); - /*What is the delta thickness forcing the sea-level rise core: cumulated over time, hence the +=:*/ + /*What is the delta thickness forcing the sea-level change core: cumulated over time, hence the +=:*/ for(int i=0;iAddBasalInput(SurfaceEnum,newsurface,P1Enum); element->AddBasalInput(BaseEnum,newbase,P1Enum); - element->AddBasalInput(SealevelriseCumDeltathicknessEnum,cumdeltathickness,P1Enum); + element->AddBasalInput(SealevelchangeCumDeltathicknessEnum,cumdeltathickness,P1Enum); element->AddBasalInput(SurfaceloadIceThicknessChangeEnum,deltathickness,P1Enum); Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 25947) +++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 25947) @@ -0,0 +1,532 @@ +#include "./SealevelchangeAnalysis.h" +#include +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../classes/Inputs/TransientInput.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +/*Model processing*/ +void SealevelchangeAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + /*No constraints*/ +}/*}}}*/ +void SealevelchangeAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + /*No loads*/ +}/*}}}*/ +void SealevelchangeAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/ + ::CreateNodes(nodes,iomodel,SealevelchangeAnalysisEnum,P1Enum); +}/*}}}*/ +int SealevelchangeAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ + return 1; +}/*}}}*/ +void SealevelchangeAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + int geodetic=0; + int dslmodel=0; + int externalnature=0; + int isexternal=0; + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(inputs,i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + /*Create inputs: */ + iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum); + iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); + iomodel->FetchData(&geodetic,"md.solidearth.settings.computesealevelchange"); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.icethicknesschange",SurfaceloadIceThicknessChangeEnum); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.waterheightchange",SurfaceloadWaterHeightChangeEnum); + + /*dynamic sea level: */ + iomodel->FetchData(&dslmodel,"md.dsl.model"); + if (dslmodel==1){ /*standard dsl model:{{{*/ + + /*deal with global mean steric rate: */ + IssmDouble* str=NULL; + IssmDouble* times = NULL; + int M,N; + + /*fetch str vector:*/ + iomodel->FetchData(&str,&M,&N,"md.dsl.global_average_thermosteric_sea_level_change"); _assert_(M==2); + + //recover time vector: + times=xNew(N); + for(int t=0;tSetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,times,N); + TransientInput* transientinput = inputs->GetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum); + + + for(Object* & object : elements->objects){ + Element* element=xDynamicCast(object); + + for(int t=0;tObjectEnum()){ + case TriaEnum: transientinput->AddTriaTimeInput( t,1,&element->lid,&str[t],P0Enum); break; + case PentaEnum: transientinput->AddPentaTimeInput(t,1,&element->lid,&str[t],P0Enum); break; + default: _error_("Not implemented yet"); + } + } + } + + /*cleanup:*/ + xDelete(times); + iomodel->DeleteData(str,"md.dsl.global_average_thermosteric_sea_level_change"); + + /*deal with dynamic sea level fields: */ + iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_surface_height_change_above_geoid", DslSeaSurfaceHeightChangeAboveGeoidEnum); + iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_water_pressure_change_at_sea_floor", DslSeaWaterPressureChangeAtSeaFloorEnum); + + } /*}}}*/ + else if (dslmodel==2){ /*multi-model ensemble dsl model:{{{*/ + + /*variables:*/ + int nummodels; + IssmDouble** pstr=NULL; + IssmDouble* str=NULL; + IssmDouble* times = NULL; + int* pM = NULL; + int* pN = NULL; + int M,N; + + /*deal with dsl.global_average_thermosteric_sea_level_change {{{*/ + iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.global_average_thermosteric_sea_level_change"); + + /*go through the mat array and create a dataset of transient inputs:*/ + for (int i=0;i(N); + for(int t=0;tSetDatasetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,i, times,N); + + for(Object* & object : elements->objects){ + Element* element=xDynamicCast(object); + + for(int t=0;tObjectEnum()){ + case TriaEnum: transientinput->AddTriaTimeInput( t,1,&element->lid,&str[t],P0Enum); break; + case PentaEnum: transientinput->AddPentaTimeInput(t,1,&element->lid,&str[t],P0Enum); break; + default: _error_("Not implemented yet"); + } + } + } + xDelete(times); + } + /*Delete data:*/ + for(int i=0;i(str); + } + xDelete(pstr); + xDelete(pM); + xDelete(pN); + /*}}}*/ + iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_surface_height_change_above_geoid",DslSeaSurfaceHeightChangeAboveGeoidEnum); + iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_water_pressure_change_at_sea_floor",DslSeaWaterPressureChangeAtSeaFloorEnum); + + } /*}}}*/ + else _error_("Dsl model " << dslmodel << " not implemented yet!"); + + /*external solidearthsolution: solid-Earth model*/ + iomodel->FetchData(&isexternal,"md.solidearth.isexternal"); + if(isexternal){ + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.external.displacementeast",SolidearthExternalDisplacementEastRateEnum); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.external.displacementnorth",SolidearthExternalDisplacementNorthRateEnum); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.external.displacementup",SolidearthExternalDisplacementUpRateEnum); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.external.geoid",SolidearthExternalGeoidRateEnum); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.external.barystaticsealevel",SolidearthExternalBarystaticSeaLevelRateEnum); + } + + /*Initialize cumdeltalthickness, bedrockeast and bedrocknorth to 0 at time step 0:*/ + InputUpdateFromConstantx(inputs,elements,0.,BedEastEnum); + InputUpdateFromConstantx(inputs,elements,0.,BedNorthEnum); + InputUpdateFromConstantx(inputs,elements,0.,SealevelchangeCumDeltathicknessEnum); + + +}/*}}}*/ +void SealevelchangeAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + int nl; + IssmDouble* love_h=NULL; + IssmDouble* love_k=NULL; + IssmDouble* love_l=NULL; + IssmDouble* love_th=NULL; + IssmDouble* love_tk=NULL; + IssmDouble* love_tl=NULL; + int dslmodel=0; + int externalnature=0; + + IssmDouble* G_rigid = NULL; + IssmDouble* G_rigid_local = NULL; + IssmDouble* G_elastic = NULL; + IssmDouble* G_elastic_local = NULL; + IssmDouble* U_elastic = NULL; + IssmDouble* U_elastic_local = NULL; + IssmDouble* H_elastic = NULL; + IssmDouble* H_elastic_local = NULL; + int M,m,lower_row,upper_row; + IssmDouble degacc=.01; + IssmDouble planetradius=0; + IssmDouble planetarea=0; + bool rigid=false; + bool elastic=false; + bool rotation=false; + + int numoutputs; + char** requestedoutputs = NULL; + + /*transition vectors: */ + IssmDouble **transitions = NULL; + int *transitions_M = NULL; + int *transitions_N = NULL; + int ntransitions; + IssmDouble* partitionice=NULL; + IssmDouble* partitionhydro=NULL; + IssmDouble* bslrice_partition=NULL; + IssmDouble* bslrhydro_partition=NULL; + int npartice,nparthydro,nel; + + + /*some constant parameters: */ + parameters->AddObject(iomodel->CopyConstantObject("md.dsl.model",DslModelEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.isexternal",SolidearthIsExternalEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.external.nature",SolidearthExternalNatureEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.runfrequency",SolidearthSettingsRunFrequencyEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.reltol",SolidearthSettingsReltolEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.abstol",SolidearthSettingsAbstolEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.maxiter",SolidearthSettingsMaxiterEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rigid",SolidearthSettingsRigidEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.horiz",SolidearthSettingsHorizEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.elastic",SolidearthSettingsElasticEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rotation",SolidearthSettingsRotationEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.equatorialmoi",RotationalEquatorialMoiEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.polarmoi",RotationalPolarMoiEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.angularvelocity",RotationalAngularVelocityEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.ocean_area_scaling",SolidearthSettingsOceanAreaScalingEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.computesealevelchange",SolidearthSettingsComputesealevelchangeEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.isgrd",SolidearthSettingsGRDEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.planetradius",SolidearthPlanetRadiusEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.glfraction",SolidearthSettingsGlfractionEnum)); + parameters->AddObject(new DoubleParam(CumBslrEnum,0.0)); + parameters->AddObject(new DoubleParam(CumBslrIceEnum,0.0)); + parameters->AddObject(new DoubleParam(CumBslrHydroEnum,0.0)); + parameters->AddObject(new DoubleParam(CumGmtslrEnum,0.0)); + + /*compute planet area and plug into parameters:*/ + iomodel->FetchData(&planetradius,"md.solidearth.planetradius"); + planetarea=4*PI*planetradius*planetradius; + parameters->AddObject(new DoubleParam(SolidearthPlanetAreaEnum,planetarea)); + + /*Deal with partition of the barystatic contribution:*/ + iomodel->FetchData(&npartice,"md.solidearth.npartice"); + parameters->AddObject(new IntParam(SolidearthNpartIceEnum,npartice)); + if(npartice){ + iomodel->FetchData(&partitionice,&nel,NULL,"md.solidearth.partitionice"); + parameters->AddObject(new DoubleMatParam(SolidearthPartitionIceEnum,partitionice,nel,1)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.npartice",SolidearthNpartIceEnum)); + bslrice_partition=xNewZeroInit(npartice); + parameters->AddObject(new DoubleMatParam(CumBslrIcePartitionEnum,bslrice_partition,npartice,1)); + xDelete(partitionice); + } + iomodel->FetchData(&nparthydro,"md.solidearth.nparthydro"); + parameters->AddObject(new IntParam(SolidearthNpartHydroEnum,nparthydro)); + if(nparthydro){ + iomodel->FetchData(&partitionhydro,&nel,NULL,"md.solidearth.partitionhydro"); + parameters->AddObject(new DoubleMatParam(SolidearthPartitionHydroEnum,partitionhydro,nel,1)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.nparthydro",SolidearthNpartHydroEnum)); + bslrhydro_partition=xNewZeroInit(nparthydro); + parameters->AddObject(new DoubleMatParam(CumBslrHydroPartitionEnum,bslrhydro_partition,nparthydro,1)); + xDelete(partitionhydro); + } + + /*Deal with dsl multi-model ensembles: {{{*/ + iomodel->FetchData(&dslmodel,"md.dsl.model"); + parameters->AddObject(iomodel->CopyConstantObject("md.dsl.compute_fingerprints",DslComputeFingerprintsEnum)); + if(dslmodel==2){ + IssmDouble modelid; + int nummodels; + + /*create double param, not int param, because Dakota will be updating it as + * a double potentially: */ + iomodel->FetchData(&modelid,"md.dsl.modelid"); + parameters->AddObject(new DoubleParam(DslModelidEnum,modelid)); + parameters->AddObject(iomodel->CopyConstantObject("md.dsl.nummodels",DslNummodelsEnum)); + iomodel->FetchData(&nummodels,"md.dsl.nummodels"); + + /*quick checks: */ + if(nummodels<=0)_error_("dslmme object in md.dsl field should contain at least 1 ensemble model!"); + if(modelid<=0 || modelid>nummodels)_error_("modelid field in dslmme object of md.dsl field should be between 1 and the number of ensemble runs!"); + } /*}}}*/ + /*Deal with external multi-model ensembles: {{{*/ + iomodel->FetchData(&externalnature,"md.solidearth.external.nature"); + if(externalnature>=3){ + IssmDouble modelid; + int nummodels; + + /*create double param, not int param, because Dakota will be updating it as a double potentially: */ + iomodel->FetchData(&modelid,"md.solidearth.external.modelid"); + parameters->AddObject(new DoubleParam(SolidearthExternalModelidEnum,modelid)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.external.nummodels",SolidearthExternalNummodelsEnum)); + iomodel->FetchData(&nummodels,"md.solidearth.external.nummodels"); + + /*quick checks: */ + if(nummodels<=0)_error_("mme solidearth solution object in md.solidearth.external field should contain at least 1 ensemble model!"); + if(modelid<=0 || modelid>nummodels)_error_("modelid field in md.solidearth.external field should be between 1 and the number of ensemble runs!"); + } /*}}}*/ + + /*Deal with elasticity {{{*/ + iomodel->FetchData(&rigid,"md.solidearth.settings.rigid"); + iomodel->FetchData(&elastic,"md.solidearth.settings.elastic"); + iomodel->FetchData(&rotation,"md.solidearth.settings.rotation"); + + if(elastic | rigid){ + /*compute green functions for a range of angles*/ + iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); + M=reCast(180./degacc+1.); + } + + /*love numbers: */ + if(elastic){ + iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h"); + iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k"); + iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l"); + iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th"); + iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk"); + iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl"); + + parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1)); + parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1)); + parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1)); + parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1)); + parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1)); + parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1)); + + // AD performance is sensitive to calls to ensurecontiguous. + // // Providing "t" will cause ensurecontiguous to be called. + #ifdef _HAVE_AD_ + G_elastic=xNew(M,"t"); + U_elastic=xNew(M,"t"); + H_elastic=xNew(M,"t"); + #else + G_elastic=xNew(M); + U_elastic=xNew(M); + H_elastic=xNew(M); + #endif + } + if(rigid){ + #ifdef _HAVE_AD_ + G_rigid=xNew(M,"t"); + #else + G_rigid=xNew(M); + #endif + } + + if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); + + if(rigid | elastic){ + + /*compute combined legendre + love number (elastic green function:*/ + m=DetermineLocalSize(M,IssmComm::GetComm()); + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); + } + if(elastic){ + #ifdef _HAVE_AD_ + G_elastic_local=xNew(m,"t"); + U_elastic_local=xNew(m,"t"); + H_elastic_local=xNew(m,"t"); + #else + G_elastic_local=xNew(m); + U_elastic_local=xNew(m); + H_elastic_local=xNew(m); + #endif + } + if(rigid){ + #ifdef _HAVE_AD_ + G_rigid_local=xNew(m,"t"); + #else + G_rigid_local=xNew(m); + #endif + } + + if(rigid){ + for(int i=lower_row;i(i)*degacc * PI / 180.0; + G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); + } + } + if(elastic){ + for(int i=lower_row;i(i)*degacc * PI / 180.0; + + G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row]; + U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row]; + H_elastic_local[i-lower_row]= 0; + IssmDouble Pn = 0.; + IssmDouble Pn1 = 0.; + IssmDouble Pn2 = 0.; + IssmDouble Pn_p = 0.; + IssmDouble Pn_p1 = 0.; + IssmDouble Pn_p2 = 0.; + + for (int n=0;n(IssmComm::GetSize()); + int* displs=xNew(IssmComm::GetSize()); + + //recvcounts: + ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); + + /*displs: */ + ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); + + /*All gather:*/ + ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + if(elastic){ + ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + } + + /*free resources: */ + xDelete(recvcounts); + xDelete(displs); + + /*Avoid singularity at 0: */ + G_rigid[0]=G_rigid[1]; + if(elastic){ + G_elastic[0]=G_elastic[1]; + U_elastic[0]=U_elastic[1]; + H_elastic[0]=H_elastic[1]; + } + + + /*Save our precomputed tables into parameters*/ + parameters->AddObject(new DoubleVecParam(SealevelchangeGRigidEnum,G_rigid,M)); + if(elastic){ + parameters->AddObject(new DoubleVecParam(SealevelchangeGElasticEnum,G_elastic,M)); + parameters->AddObject(new DoubleVecParam(SealevelchangeUElasticEnum,U_elastic,M)); + parameters->AddObject(new DoubleVecParam(SealevelchangeHElasticEnum,H_elastic,M)); + } + + /*free resources: */ + xDelete(G_rigid); + xDelete(G_rigid_local); + if(elastic){ + xDelete(love_h); + xDelete(love_k); + xDelete(love_l); + xDelete(love_th); + xDelete(love_tk); + xDelete(love_tl); + xDelete(G_elastic); + xDelete(G_elastic_local); + xDelete(U_elastic); + xDelete(U_elastic_local); + xDelete(H_elastic); + xDelete(H_elastic_local); + } + } /*}}}*/ + + /*Indicate we have not yet run the Geometry Core module: */ + parameters->AddObject(new BoolParam(SealevelchangeGeometryDoneEnum,false)); + /*}}}*/ + + /*Transitions:{{{ */ + iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.solidearth.transitions"); + if(transitions){ + parameters->AddObject(new DoubleMatArrayParam(SealevelchangeTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N)); + + for(int i=0;i(transition); + } + xDelete(transitions); + xDelete(transitions_M); + xDelete(transitions_N); + } /*}}}*/ + /*Requested outputs {{{*/ + iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.solidearth.requested_outputs"); + if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelchangeRequestedOutputsEnum,requestedoutputs,numoutputs)); + iomodel->DeleteData(&requestedoutputs,numoutputs,"md.solidearth.requested_outputs"); + /*}}}*/ + +}/*}}}*/ + +/*Finite Element Analysis*/ +void SealevelchangeAnalysis::Core(FemModel* femmodel){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +ElementVector* SealevelchangeAnalysis::CreateDVector(Element* element){/*{{{*/ + /*Default, return NULL*/ + return NULL; +}/*}}}*/ +ElementMatrix* SealevelchangeAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ +_error_("Not implemented"); +}/*}}}*/ +ElementMatrix* SealevelchangeAnalysis::CreateKMatrix(Element* element){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +ElementVector* SealevelchangeAnalysis::CreatePVector(Element* element){/*{{{*/ +_error_("not implemented yet"); +}/*}}}*/ +void SealevelchangeAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void SealevelchangeAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/ + _error_("Not implemented yet"); +}/*}}}*/ +void SealevelchangeAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + _error_("not implemeneted yet!"); + +}/*}}}*/ +void SealevelchangeAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ + /*Default, do nothing*/ + return; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.h (revision 25947) +++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.h (revision 25947) @@ -0,0 +1,33 @@ +/*! \file SealevelchangeAnalysis.h + * \brief: header file for generic external result object + */ + +#ifndef _SealevelchangeAnalysis_ +#define _SealevelchangeAnalysis_ + +/*Headers*/ +#include "./Analysis.h" + +class SealevelchangeAnalysis: public Analysis{ + + public: + /*Model processing*/ + void CreateConstraints(Constraints* constraints,IoModel* iomodel); + void CreateLoads(Loads* loads, IoModel* iomodel); + void CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr=false); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); + + /*Finite element Analysis*/ + void Core(FemModel* femmodel); + ElementVector* CreateDVector(Element* element); + ElementMatrix* CreateJacobianMatrix(Element* element); + ElementMatrix* CreateKMatrix(Element* element); + ElementVector* CreatePVector(Element* element); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_interp,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); +}; +#endif Index: sm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 25946) +++ (revision ) @@ -1,504 +1,0 @@ -#include "./SealevelriseAnalysis.h" -#include -#include "../toolkits/toolkits.h" -#include "../classes/classes.h" -#include "../classes/Inputs/TransientInput.h" -#include "../shared/shared.h" -#include "../modules/modules.h" - -/*Model processing*/ -void SealevelriseAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ - /*No constraints*/ -}/*}}}*/ -void SealevelriseAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ - /*No loads*/ -}/*}}}*/ -void SealevelriseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/ - ::CreateNodes(nodes,iomodel,SealevelriseAnalysisEnum,P1Enum); -}/*}}}*/ -int SealevelriseAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ - return 1; -}/*}}}*/ -void SealevelriseAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ - - int geodetic=0; - int dslmodel=0; - - /*Update elements: */ - int counter=0; - for(int i=0;inumberofelements;i++){ - if(iomodel->my_elements[i]){ - Element* element=(Element*)elements->GetObjectByOffset(counter); - element->Update(inputs,i,iomodel,analysis_counter,analysis_type,P1Enum); - counter++; - } - } - - /*Create inputs: */ - iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum); - iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); - iomodel->FetchData(&geodetic,"md.solidearth.settings.computesealevelchange"); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.icethicknesschange",SurfaceloadIceThicknessChangeEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); - iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.waterheightchange",SurfaceloadWaterHeightChangeEnum); - - /*dynamic sea level: */ - iomodel->FetchData(&dslmodel,"md.dsl.model"); - if (dslmodel==1){ /*standard dsl model:{{{*/ - - /*deal with global mean steric rate: */ - IssmDouble* str=NULL; - IssmDouble* times = NULL; - int M,N; - - /*fetch str vector:*/ - iomodel->FetchData(&str,&M,&N,"md.dsl.global_average_thermosteric_sea_level_change"); _assert_(M==2); - - //recover time vector: - times=xNew(N); - for(int t=0;tSetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,times,N); - TransientInput* transientinput = inputs->GetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum); - - - for(Object* & object : elements->objects){ - Element* element=xDynamicCast(object); - - for(int t=0;tObjectEnum()){ - case TriaEnum: transientinput->AddTriaTimeInput( t,1,&element->lid,&str[t],P0Enum); break; - case PentaEnum: transientinput->AddPentaTimeInput(t,1,&element->lid,&str[t],P0Enum); break; - default: _error_("Not implemented yet"); - } - } - } - - /*cleanup:*/ - xDelete(times); - iomodel->DeleteData(str,"md.dsl.global_average_thermosteric_sea_level_change"); - - /*deal with dynamic sea level fields: */ - iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_surface_height_change_above_geoid", DslSeaSurfaceHeightChangeAboveGeoidEnum); - iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_water_pressure_change_at_sea_floor", DslSeaWaterPressureChangeAtSeaFloorEnum); - - } /*}}}*/ - else if (dslmodel==2){ /*multi-model ensemble dsl model:{{{*/ - - /*variables:*/ - int nummodels; - IssmDouble** pstr=NULL; - IssmDouble* str=NULL; - IssmDouble* times = NULL; - int* pM = NULL; - int* pN = NULL; - int M,N; - - /*deal with dsl.global_average_thermosteric_sea_level_change {{{*/ - iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.global_average_thermosteric_sea_level_change"); - - /*go through the mat array and create a dataset of transient inputs:*/ - for (int i=0;i(N); - for(int t=0;tSetDatasetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,i, times,N); - - for(Object* & object : elements->objects){ - Element* element=xDynamicCast(object); - - for(int t=0;tObjectEnum()){ - case TriaEnum: transientinput->AddTriaTimeInput( t,1,&element->lid,&str[t],P0Enum); break; - case PentaEnum: transientinput->AddPentaTimeInput(t,1,&element->lid,&str[t],P0Enum); break; - default: _error_("Not implemented yet"); - } - } - } - xDelete(times); - } - /*Delete data:*/ - for(int i=0;i(str); - } - xDelete(pstr); - xDelete(pM); - xDelete(pN); - /*}}}*/ - iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_surface_height_change_above_geoid",DslSeaSurfaceHeightChangeAboveGeoidEnum); - iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_water_pressure_change_at_sea_floor",DslSeaWaterPressureChangeAtSeaFloorEnum); - - } /*}}}*/ - else _error_("Dsl model " << dslmodel << " not implemented yet!"); - - - - /*Initialize cumdeltalthickness and sealevel rise rate input*/ - InputUpdateFromConstantx(inputs,elements,0.,SealevelriseCumDeltathicknessEnum); - InputUpdateFromConstantx(inputs,elements,0.,SealevelNEsaRateEnum); - InputUpdateFromConstantx(inputs,elements,0.,SealevelUEsaRateEnum); - InputUpdateFromConstantx(inputs,elements,0.,SealevelRSLRateEnum); - InputUpdateFromConstantx(inputs,elements,0.,SealevelEustaticMaskEnum); - InputUpdateFromConstantx(inputs,elements,0.,SealevelEustaticOceanMaskEnum); - -}/*}}}*/ -void SealevelriseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ - - int nl; - IssmDouble* love_h=NULL; - IssmDouble* love_k=NULL; - IssmDouble* love_l=NULL; - IssmDouble* love_th=NULL; - IssmDouble* love_tk=NULL; - IssmDouble* love_tl=NULL; - int dslmodel=0; - - IssmDouble* G_rigid = NULL; - IssmDouble* G_rigid_local = NULL; - IssmDouble* G_elastic = NULL; - IssmDouble* G_elastic_local = NULL; - IssmDouble* U_elastic = NULL; - IssmDouble* U_elastic_local = NULL; - IssmDouble* H_elastic = NULL; - IssmDouble* H_elastic_local = NULL; - int M,m,lower_row,upper_row; - IssmDouble degacc=.01; - IssmDouble planetradius=0; - IssmDouble planetarea=0; - bool rigid=false; - bool elastic=false; - bool rotation=false; - - int numoutputs; - char** requestedoutputs = NULL; - - /*transition vectors: */ - IssmDouble **transitions = NULL; - int *transitions_M = NULL; - int *transitions_N = NULL; - int ntransitions; - IssmDouble* partitionice=NULL; - IssmDouble* partitionhydro=NULL; - IssmDouble* bslrice_partition=NULL; - IssmDouble* bslrhydro_partition=NULL; - int npartice,nparthydro,nel; - - - /*some constant parameters: */ - parameters->AddObject(iomodel->CopyConstantObject("md.dsl.model",DslModelEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.runfrequency",SolidearthSettingsRunFrequencyEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.reltol",SolidearthSettingsReltolEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.abstol",SolidearthSettingsAbstolEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.maxiter",SolidearthSettingsMaxiterEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rigid",SolidearthSettingsRigidEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.horiz",SolidearthSettingsHorizEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.elastic",SolidearthSettingsElasticEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rotation",SolidearthSettingsRotationEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.equatorialmoi",RotationalEquatorialMoiEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.polarmoi",RotationalPolarMoiEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.angularvelocity",RotationalAngularVelocityEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.ocean_area_scaling",SolidearthSettingsOceanAreaScalingEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.computesealevelchange",SolidearthSettingsComputesealevelchangeEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.isgrd",SolidearthSettingsGRDEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.planetradius",SolidearthPlanetRadiusEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.glfraction",SolidearthSettingsGlfractionEnum)); - parameters->AddObject(new DoubleParam(CumBslrEnum,0.0)); - parameters->AddObject(new DoubleParam(CumBslrIceEnum,0.0)); - parameters->AddObject(new DoubleParam(CumBslrHydroEnum,0.0)); - parameters->AddObject(new DoubleParam(CumGmtslrEnum,0.0)); - - /*compute planet area and plug into parameters:*/ - iomodel->FetchData(&planetradius,"md.solidearth.planetradius"); - planetarea=4*PI*planetradius*planetradius; - parameters->AddObject(new DoubleParam(SolidearthPlanetAreaEnum,planetarea)); - - /*Deal with partition of the barystatic contribution:*/ - iomodel->FetchData(&npartice,"md.solidearth.npartice"); - parameters->AddObject(new IntParam(SolidearthNpartIceEnum,npartice)); - if(npartice){ - iomodel->FetchData(&partitionice,&nel,NULL,"md.solidearth.partitionice"); - parameters->AddObject(new DoubleMatParam(SolidearthPartitionIceEnum,partitionice,nel,1)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.npartice",SolidearthNpartIceEnum)); - bslrice_partition=xNewZeroInit(npartice); - parameters->AddObject(new DoubleMatParam(CumBslrIcePartitionEnum,bslrice_partition,npartice,1)); - xDelete(partitionice); - } - iomodel->FetchData(&nparthydro,"md.solidearth.nparthydro"); - parameters->AddObject(new IntParam(SolidearthNpartHydroEnum,nparthydro)); - if(nparthydro){ - iomodel->FetchData(&partitionhydro,&nel,NULL,"md.solidearth.partitionhydro"); - parameters->AddObject(new DoubleMatParam(SolidearthPartitionHydroEnum,partitionhydro,nel,1)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.nparthydro",SolidearthNpartHydroEnum)); - bslrhydro_partition=xNewZeroInit(nparthydro); - parameters->AddObject(new DoubleMatParam(CumBslrHydroPartitionEnum,bslrhydro_partition,nparthydro,1)); - xDelete(partitionhydro); - } - - /*Deal with dsl multi-model ensembles: {{{*/ - iomodel->FetchData(&dslmodel,"md.dsl.model"); - parameters->AddObject(iomodel->CopyConstantObject("md.dsl.compute_fingerprints",DslComputeFingerprintsEnum)); - if(dslmodel==2){ - IssmDouble modelid; - int nummodels; - - /*create double param, not int param, because Dakota will be updating it as - * a double potentially: */ - iomodel->FetchData(&modelid,"md.dsl.modelid"); - parameters->AddObject(new DoubleParam(DslModelidEnum,modelid)); - parameters->AddObject(iomodel->CopyConstantObject("md.dsl.nummodels",DslNummodelsEnum)); - iomodel->FetchData(&nummodels,"md.dsl.nummodels"); - - /*quick checks: */ - if(nummodels<=0)_error_("dslmme object in md.dsl field should contain at least 1 ensemble model!"); - if(modelid<=0 || modelid>nummodels)_error_("modelid field in dslmme object of md.dsl field should be between 1 and the number of ensemble runs!"); - } /*}}}*/ - /*Deal with elasticity {{{*/ - iomodel->FetchData(&rigid,"md.solidearth.settings.rigid"); - iomodel->FetchData(&elastic,"md.solidearth.settings.elastic"); - iomodel->FetchData(&rotation,"md.solidearth.settings.rotation"); - - if(elastic | rigid){ - /*compute green functions for a range of angles*/ - iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); - M=reCast(180./degacc+1.); - } - - /*love numbers: */ - if(elastic){ - iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h"); - iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k"); - iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l"); - iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th"); - iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk"); - iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl"); - - parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1)); - parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1)); - parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1)); - parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1)); - parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1)); - parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1)); - - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. - #ifdef _HAVE_AD_ - G_elastic=xNew(M,"t"); - U_elastic=xNew(M,"t"); - H_elastic=xNew(M,"t"); - #else - G_elastic=xNew(M); - U_elastic=xNew(M); - H_elastic=xNew(M); - #endif - } - if(rigid){ - #ifdef _HAVE_AD_ - G_rigid=xNew(M,"t"); - #else - G_rigid=xNew(M); - #endif - } - - if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); - - if(rigid | elastic){ - - /*compute combined legendre + love number (elastic green function:*/ - m=DetermineLocalSize(M,IssmComm::GetComm()); - GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); - } - if(elastic){ - #ifdef _HAVE_AD_ - G_elastic_local=xNew(m,"t"); - U_elastic_local=xNew(m,"t"); - H_elastic_local=xNew(m,"t"); - #else - G_elastic_local=xNew(m); - U_elastic_local=xNew(m); - H_elastic_local=xNew(m); - #endif - } - if(rigid){ - #ifdef _HAVE_AD_ - G_rigid_local=xNew(m,"t"); - #else - G_rigid_local=xNew(m); - #endif - } - - if(rigid){ - for(int i=lower_row;i(i)*degacc * PI / 180.0; - G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); - } - } - if(elastic){ - for(int i=lower_row;i(i)*degacc * PI / 180.0; - - G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row]; - U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row]; - H_elastic_local[i-lower_row]= 0; - IssmDouble Pn = 0.; - IssmDouble Pn1 = 0.; - IssmDouble Pn2 = 0.; - IssmDouble Pn_p = 0.; - IssmDouble Pn_p1 = 0.; - IssmDouble Pn_p2 = 0.; - - for (int n=0;n(IssmComm::GetSize()); - int* displs=xNew(IssmComm::GetSize()); - - //recvcounts: - ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); - - /*displs: */ - ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); - - /*All gather:*/ - ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - if(elastic){ - ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - } - - /*free resources: */ - xDelete(recvcounts); - xDelete(displs); - - /*Avoid singularity at 0: */ - G_rigid[0]=G_rigid[1]; - if(elastic){ - G_elastic[0]=G_elastic[1]; - U_elastic[0]=U_elastic[1]; - H_elastic[0]=H_elastic[1]; - } - - - /*Save our precomputed tables into parameters*/ - parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); - if(elastic){ - parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M)); - parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M)); - parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M)); - } - - /*free resources: */ - xDelete(G_rigid); - xDelete(G_rigid_local); - if(elastic){ - xDelete(love_h); - xDelete(love_k); - xDelete(love_l); - xDelete(love_th); - xDelete(love_tk); - xDelete(love_tl); - xDelete(G_elastic); - xDelete(G_elastic_local); - xDelete(U_elastic); - xDelete(U_elastic_local); - xDelete(H_elastic); - xDelete(H_elastic_local); - } - } /*}}}*/ - - /*Indicate we have not yet run the Geometry Core module: */ - parameters->AddObject(new BoolParam(SealevelriseGeometryDoneEnum,false)); - /*}}}*/ - - /*Transitions:{{{ */ - iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.solidearth.transitions"); - if(transitions){ - parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N)); - - for(int i=0;i(transition); - } - xDelete(transitions); - xDelete(transitions_M); - xDelete(transitions_N); - } /*}}}*/ - /*Requested outputs {{{*/ - iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.solidearth.requested_outputs"); - if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs)); - iomodel->DeleteData(&requestedoutputs,numoutputs,"md.solidearth.requested_outputs"); - /*}}}*/ - -}/*}}}*/ - -/*Finite Element Analysis*/ -void SealevelriseAnalysis::Core(FemModel* femmodel){/*{{{*/ - _error_("not implemented"); -}/*}}}*/ -ElementVector* SealevelriseAnalysis::CreateDVector(Element* element){/*{{{*/ - /*Default, return NULL*/ - return NULL; -}/*}}}*/ -ElementMatrix* SealevelriseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ -_error_("Not implemented"); -}/*}}}*/ -ElementMatrix* SealevelriseAnalysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ -ElementVector* SealevelriseAnalysis::CreatePVector(Element* element){/*{{{*/ -_error_("not implemented yet"); -}/*}}}*/ -void SealevelriseAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ -void SealevelriseAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/ - _error_("Not implemented yet"); -}/*}}}*/ -void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - _error_("not implemeneted yet!"); - -}/*}}}*/ -void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ - /*Default, do nothing*/ - return; -}/*}}}*/ Index: sm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.h (revision 25946) +++ (revision ) @@ -1,33 +1,0 @@ -/*! \file SealevelriseAnalysis.h - * \brief: header file for generic external result object - */ - -#ifndef _SealevelriseAnalysis_ -#define _SealevelriseAnalysis_ - -/*Headers*/ -#include "./Analysis.h" - -class SealevelriseAnalysis: public Analysis{ - - public: - /*Model processing*/ - void CreateConstraints(Constraints* constraints,IoModel* iomodel); - void CreateLoads(Loads* loads, IoModel* iomodel); - void CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr=false); - int DofsPerNode(int** doflist,int domaintype,int approximation); - void UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type); - void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); - - /*Finite element Analysis*/ - void Core(FemModel* femmodel); - ElementVector* CreateDVector(Element* element); - ElementMatrix* CreateJacobianMatrix(Element* element); - ElementMatrix* CreateKMatrix(Element* element); - ElementVector* CreatePVector(Element* element); - void GetSolutionFromInputs(Vector* solution,Element* element); - void GradientJ(Vector* gradient,Element* element,int control_type,int control_interp,int control_index); - void InputUpdateFromSolution(IssmDouble* solution,Element* element); - void UpdateConstraints(FemModel* femmodel); -}; -#endif Index: /issm/trunk-jpl/src/c/analyses/analyses.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/analyses.h (revision 25946) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 25947) @@ -36,5 +36,5 @@ #include "./MasstransportAnalysis.h" #include "./SmbAnalysis.h" -#include "./SealevelriseAnalysis.h" +#include "./SealevelchangeAnalysis.h" #include "./MeltingAnalysis.h" #include "./SmoothAnalysis.h" Index: /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 25946) +++ /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 25947) @@ -95,5 +95,5 @@ /*Hack:*/ - for(int i=0;inummodels;i++) if(femmodel_init->analysis_type_list[i]==SealevelriseAnalysisEnum) sealevelrise_core_geometry(femmodel_init); + for(int i=0;inummodels;i++) if(femmodel_init->analysis_type_list[i]==SealevelchangeAnalysisEnum) sealevelchange_geometry(femmodel_init); /*Make a copy of femmodel, so we start this new evaluation run for this specific sample with a brand Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 25946) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 25947) @@ -377,11 +377,11 @@ virtual IssmDouble GetAreaSpherical(void)=0; virtual IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks)=0; - virtual void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks)=0; - virtual IssmDouble SealevelriseEustaticIce(IssmDouble* Sgi, SealevelMasks* masks,Vector* barystatic_contribution,IssmDouble* partition,IssmDouble oceanarea)=0; - virtual IssmDouble SealevelriseEustaticHydro(IssmDouble* Sgi, SealevelMasks* masks,Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea)=0; - virtual void SealevelriseEustaticBottomPressure(IssmDouble* Sgi, SealevelMasks* masks)=0; - virtual void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz)=0; - virtual void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* mask)=0; - virtual void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks)=0; + virtual void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks)=0; + virtual IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks,Vector* barystatic_contribution,IssmDouble* partition,IssmDouble oceanarea)=0; + virtual IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks,Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea)=0; + virtual void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi, SealevelMasks* masks)=0; + virtual void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz)=0; + virtual void SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* mask)=0; + virtual void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 25946) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 25947) @@ -216,11 +216,11 @@ IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; - IssmDouble SealevelriseEustaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");}; - IssmDouble SealevelriseEustaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");}; - void SealevelriseEustaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; + IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");}; + IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");}; + void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; + void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 25946) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 25947) @@ -170,12 +170,12 @@ #endif #ifdef _HAVE_SEALEVELRISE_ - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; - IssmDouble SealevelriseEustaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");}; - IssmDouble SealevelriseEustaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");}; - void SealevelriseEustaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; + IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");}; + IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");}; + void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; + void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 25946) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 25947) @@ -177,11 +177,11 @@ #ifdef _HAVE_SEALEVELRISE_ void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; - IssmDouble SealevelriseEustaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");}; - IssmDouble SealevelriseEustaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");}; - void SealevelriseEustaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; + IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");}; + IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");}; + void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; + void DeformationFromSurfaceLoads(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25946) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25947) @@ -5387,5 +5387,5 @@ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); - /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ + /*Get area of element: precomputed in the sealevelchange_geometry:*/ area=GetAreaSpherical(); @@ -5522,5 +5522,5 @@ } /*}}}*/ -void Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/ +void Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/ /*early return if we are not on an ice cap OR ocean:*/ if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){ @@ -5586,5 +5586,5 @@ rho_water=FindParam(MaterialsRhoSeawaterEnum); - /*From Sg_old, recover water sea level rise:*/ + /*From Sg_old, recover water sea level change:*/ S=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; @@ -5606,5 +5606,5 @@ /*Compute ice thickness change: */ Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessChangeEnum); - if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); + if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!"); deltathickness_input->GetInputAverage(&I); @@ -5632,5 +5632,5 @@ } /*}}}*/ -void Tria::SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){ /*{{{*/ +void Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){ /*{{{*/ /*diverse:*/ int gsize; @@ -5693,5 +5693,5 @@ /*recover precomputed green function kernels:*/ - DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); + DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGRigidEnum)); _assert_(parameter); parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M); @@ -5701,11 +5701,11 @@ if(computeelastic){ - parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); + parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M); - parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter); + parameter = static_cast(this->parameters->FindParamObject(SealevelchangeHElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M); - parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter); + parameter = static_cast(this->parameters->FindParamObject(SealevelchangeUElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); @@ -5775,11 +5775,11 @@ /*Add in inputs:*/ - this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); - this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); + this->inputs->SetArrayInput(SealevelchangeIndicesEnum,this->lid,indices,gsize); + this->inputs->SetArrayInput(SealevelchangeGEnum,this->lid,G,gsize); if(computeelastic){ - this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize); + this->inputs->SetArrayInput(SealevelchangeGUEnum,this->lid,GU,gsize); if(horiz){ - this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize); - this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize); + this->inputs->SetArrayInput(SealevelchangeGNEnum,this->lid,GN,gsize); + this->inputs->SetArrayInput(SealevelchangeGEEnum,this->lid,GE,gsize); } } @@ -5813,10 +5813,10 @@ } /*}}}*/ -IssmDouble Tria::SealevelriseEustaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/ +IssmDouble Tria::SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/ /*diverse:*/ int gsize; IssmDouble area; - IssmDouble phi=1.0; //WARNING: do not touch this, default is entire elemnt contributes eustatic + IssmDouble phi=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) bool notfullygrounded=false; @@ -5840,5 +5840,5 @@ if(!masks->isiceonly[this->lid]){ #ifdef _ISSM_DEBUG_ - constant=0; this->AddInput(SealevelEustaticMaskEnum,&constant,P0Enum); + constant=0; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); #endif bslrice=0; @@ -5850,5 +5850,5 @@ constant=0; #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelEustaticMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); #endif bslrice=0; @@ -5862,5 +5862,5 @@ constant=1; #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelEustaticMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); #endif @@ -5873,7 +5873,7 @@ /*retrieve precomputed G:*/ - if(computerigid)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); - - /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ + if(computerigid)this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&gsize); + + /*Get area of element: precomputed in the sealevelchange_geometry:*/ this->Element::GetInputValue(&area,AreaEnum); @@ -5886,5 +5886,5 @@ if(glfraction==0)phi=1; #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelEustaticMaskEnum,&phi,P0Enum); + this->AddInput(SealevelBarystaticMaskEnum,&phi,P0Enum); #endif } @@ -5893,5 +5893,5 @@ /*Retrieve ice thickness at vertices: */ Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessChangeEnum); - if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); + if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!"); /*/Average ice thickness over grounded area of the element only: {{{*/ @@ -5944,5 +5944,5 @@ } /*}}}*/ -IssmDouble Tria::SealevelriseEustaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/ +IssmDouble Tria::SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/ /*diverse:*/ @@ -5987,12 +5987,12 @@ /*retrieve precomputed G:*/ - if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); - - /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ + if(computeelastic)this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&gsize); + + /*Get area of element: precomputed in the sealevelchange_geometry:*/ this->Element::GetInputValue(&area,AreaEnum); /*Retrieve water height at vertices: */ Input* deltathickness_input=this->GetInput(SurfaceloadWaterHeightChangeEnum); - if (!deltathickness_input)_error_("SurfaceloadWaterHeightChangeEnum input needed to compute sea level rise!"); + if (!deltathickness_input)_error_("SurfaceloadWaterHeightChangeEnum input needed to compute sea level change!"); deltathickness_input->GetInputAverage(&W); @@ -6020,5 +6020,5 @@ } /*}}}*/ -void Tria::SealevelriseEustaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/ +void Tria::SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/ /*diverse:*/ @@ -6044,5 +6044,5 @@ /*Inform mask: */ #ifdef _ISSM_DEBUG_ - constant=1; this->AddInput(SealevelEustaticMaskEnum,&constant,P0Enum); + constant=1; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); #endif @@ -6051,12 +6051,12 @@ /*retrieve precomputed G:*/ - this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); - - /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ + this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&gsize); + + /*Get area of element: precomputed in the sealevelchange_geometry:*/ this->Element::GetInputValue(&area,AreaEnum); /*Retrieve bottom pressure change and average over the element: */ Input* bottompressure_change_input=this->GetInput(DslSeaWaterPressureChangeAtSeaFloorEnum); - if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level rise fingerprint!"); + if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level change fingerprint!"); bottompressure_change_input->GetInputAverage(&BP); @@ -6070,5 +6070,5 @@ } /*}}}*/ -void Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){ /*{{{*/ +void Tria::SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){ /*{{{*/ /*diverse:*/ @@ -6088,5 +6088,5 @@ constant=0; #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelEustaticOceanMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticOceanMaskEnum,&constant,P0Enum); #endif return; @@ -6094,5 +6094,5 @@ constant=1; #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelEustaticOceanMaskEnum,&constant,P0Enum); + this->AddInput(SealevelBarystaticOceanMaskEnum,&constant,P0Enum); #endif @@ -6101,7 +6101,7 @@ /*retrieve precomputed G:*/ - this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&dummy); _assert_(dummy==gsize); - - /*From Sg_old, recover water sea level rise:*/ + this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&dummy); _assert_(dummy==gsize); + + /*From Sg_old, recover water sea level change:*/ S=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; @@ -6114,5 +6114,5 @@ } /*}}}*/ -void Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/ +void Tria::DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/ /*diverse:*/ @@ -6147,13 +6147,13 @@ /*recover elastic Green's functions for displacement:*/ - this->inputs->GetArrayPtr(SealevelriseGUEnum,this->lid,&GU,&gsize); + this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&gsize); if(horiz){ - this->inputs->GetArrayPtr(SealevelriseGEEnum,this->lid,&GE,&gsize); - this->inputs->GetArrayPtr(SealevelriseGNEnum,this->lid,&GN,&gsize); + this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&gsize); + this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&gsize); } if(masks->isoceanin[this->lid]){ - /*From Sg, recover water sea level rise:*/ + /*From Sg, recover water sea level change:*/ S=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; @@ -6164,5 +6164,5 @@ if(bp_compute_fingerprints){ Input* bottompressure_change_input=this->GetInput(DslSeaWaterPressureChangeAtSeaFloorEnum); - if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level rise fingerprint!"); + if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level change fingerprint!"); bottompressure_change_input->GetInputAverage(&BP); @@ -6185,5 +6185,5 @@ /*Compute ice thickness change: */ Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessChangeEnum); - if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); + if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!"); deltathickness_input->GetInputAverage(&I); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 25946) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 25947) @@ -162,12 +162,12 @@ #ifdef _HAVE_SEALEVELRISE_ IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks); - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks); + void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks); void SetSealevelMasks(SealevelMasks* masks); - void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); - IssmDouble SealevelriseEustaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea); - IssmDouble SealevelriseEustaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea); - void SealevelriseEustaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks); - void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks); - void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks); + void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); + IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea); + IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea); + void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks); + void SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks); + void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks); #endif /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 25946) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 25947) @@ -833,6 +833,6 @@ break; - case SealevelriseSolutionEnum: - analyses_temp[numanalyses++]=SealevelriseAnalysisEnum; + case SealevelchangeSolutionEnum: + analyses_temp[numanalyses++]=SealevelchangeAnalysisEnum; analyses_temp[numanalyses++]=GiaAnalysisEnum; break; @@ -899,5 +899,5 @@ if(isslr){ analyses_temp[numanalyses++]=GiaAnalysisEnum; - analyses_temp[numanalyses++]=SealevelriseAnalysisEnum; + analyses_temp[numanalyses++]=SealevelchangeAnalysisEnum; } if(isesa){ @@ -4757,5 +4757,5 @@ #endif #ifdef _HAVE_SEALEVELRISE_ -void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslr,IssmDouble* pbslrice, IssmDouble* pbslrhydro, IssmDouble** pbslrice_partition,IssmDouble** pbslrhydro_partition,SealevelMasks* masks) { /*{{{*/ +void FemModel::SealevelchangeBarystatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslr,IssmDouble* pbslrice, IssmDouble* pbslrhydro, IssmDouble** pbslrice_partition,IssmDouble** pbslrhydro_partition,SealevelMasks* masks) { /*{{{*/ /*serialized vectors:*/ @@ -4780,5 +4780,5 @@ - /*Initialize temporary vector that will be used to sum eustatic components + /*Initialize temporary vector that will be used to sum barystatic components * on all local elements, prior to assembly:*/ int gsize = this->nodes->NumberOfDofs(GsetEnum); @@ -4787,5 +4787,5 @@ for(int i=0;ielements->objects){ @@ -4813,24 +4813,24 @@ - /*Call the sea level rise core for ice : */ + /*Call the barystatic sea level change core for ice : */ bslrice_cpu=0; for(Object* & object : this->elements->objects){ Element* element = xDynamicCast(object); - bslrice_cpu+=element->SealevelriseEustaticIce(RSLgi,masks, bslrice_partition,partitionice,oceanarea); - } - - /*Call the sea level rise core for hydro: */ + bslrice_cpu+=element->SealevelchangeBarystaticIce(RSLgi,masks, bslrice_partition,partitionice,oceanarea); + } + + /*Call the barystatic sea level change core for hydro: */ bslrhydro_cpu=0; for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - bslrhydro_cpu+=element->SealevelriseEustaticHydro(RSLgi,masks, bslrhydro_partition,partitionhydro,oceanarea); - } - - /*Call the sea level rise core for bottom pressures: */ + bslrhydro_cpu+=element->SealevelchangeBarystaticHydro(RSLgi,masks, bslrhydro_partition,partitionhydro,oceanarea); + } + + /*Call the barystatic sea level change core for bottom pressures: */ this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum); if(bp_compute_fingerprints){ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseEustaticBottomPressure(RSLgi,masks); + element->SealevelchangeBarystaticBottomPressure(RSLgi,masks); } } @@ -4840,5 +4840,5 @@ pRSLgi->Assemble(); - /*Sum all eustatic components from all cpus:*/ + /*Sum all barystatic components from all cpus:*/ ISSM_MPI_Reduce (&bslrice_cpu,&bslrice,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&bslrice,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); @@ -4878,5 +4878,5 @@ } /*}}}*/ -void FemModel::SealevelriseNonEustatic(Vector* pRSLgo, Vector* pRSLg_old, SealevelMasks* masks, bool verboseconvolution){/*{{{*/ +void FemModel::SealevelchangeSal(Vector* pRSLgo, Vector* pRSLg_old, SealevelMasks* masks, bool verboseconvolution){/*{{{*/ /*serialized vectors:*/ @@ -4893,5 +4893,5 @@ this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); - /*Initialize temporary vector that will be used to sum eustatic components on all local elements, prior + /*Initialize temporary vector that will be used to sum barystatic components on all local elements, prior * to assembly:*/ gsize = this->nodes->NumberOfDofs(GsetEnum); @@ -4902,9 +4902,9 @@ RSLg_old=pRSLg_old->ToMPISerial(); - /*Call the sea level rise non-eustatic core only if required: */ + /*Call the sal sea level change core only if required: */ if(computerigid){ for(Object* & object : this->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks); + element->SealevelchangeSal(RSLgo,RSLg_old,masks); } } @@ -4919,5 +4919,5 @@ } /*}}}*/ -void FemModel::SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks){/*{{{*/ +void FemModel::SealevelchangeRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks){/*{{{*/ /*serialized vectors:*/ @@ -4942,5 +4942,5 @@ for(Object* & object : this->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,masks ); + element->SealevelchangeMomentOfInertia(&moi_list[0],RSLg_old,masks ); moi_list_cpu[0] += moi_list[0]; moi_list_cpu[1] += moi_list[1]; @@ -5013,5 +5013,5 @@ } /*}}}*/ -void FemModel::SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, SealevelMasks* masks){/*{{{*/ +void FemModel::SealevelchangeDeformation(Vector* pgeoid, Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, SealevelMasks* masks){/*{{{*/ /*serialized vectors:*/ @@ -5031,5 +5031,5 @@ RSLg=pRSLg->ToMPISerial(); - /*Initialize temporary vector that will be used to sum eustatic components on all local elements, prior + /*Initialize temporary vector that will be used to sum barystatic components on all local elements, prior * to assembly:*/ gsize = this->nodes->NumberOfDofs(GsetEnum); @@ -5041,12 +5041,17 @@ indices=xNew(gsize); for (int i=0;ielements->objects){ Element* element = xDynamicCast(object); - element->SealevelriseGeodetic(Up,North,East,RSLg,masks); + element->DeformationFromSurfaceLoads(Up,North,East,RSLg,masks); } pUp->SetValues(gsize,indices,Up,ADD_VAL); pUp->Assemble(); + + + /*Add RSL to Up to find the geoid: */ + pUp->Copy(pgeoid); pgeoid->AXPY(pRSLg,1); + if (horiz){ pNorth->SetValues(gsize,indices,North,ADD_VAL); @@ -5066,5 +5071,5 @@ } /*}}}*/ -IssmDouble FemModel::SealevelriseOceanAverage(Vector* RSLg,SealevelMasks* masks, IssmDouble oceanarea) { /*{{{*/ +IssmDouble FemModel::SealevelchangeOceanAverage(Vector* RSLg,SealevelMasks* masks, IssmDouble oceanarea) { /*{{{*/ IssmDouble* RSLg_serial=NULL; Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 25946) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 25947) @@ -166,9 +166,9 @@ #endif #ifdef _HAVE_SEALEVELRISE_ - void SealevelriseEustatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslr,IssmDouble* pbslrice, IssmDouble* pbslrhydro, IssmDouble** pbslrice_partition,IssmDouble** pbslrhydro_partition, SealevelMasks* masks); - void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, SealevelMasks* masks,bool verboseconvolution); - void SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks); - void SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, SealevelMasks* masks); - IssmDouble SealevelriseOceanAverage(Vector* Sg,SealevelMasks* masks, IssmDouble oceanarea); + void SealevelchangeBarystatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslr,IssmDouble* pbslrice, IssmDouble* pbslrhydro, IssmDouble** pbslrice_partition,IssmDouble** pbslrhydro_partition, SealevelMasks* masks); + void SealevelchangeSal(Vector* pSgo, Vector* pSg_old, SealevelMasks* masks,bool verboseconvolution); + void SealevelchangeRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks); + void SealevelchangeDeformation(Vector* pgeoid,Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, SealevelMasks* masks); + IssmDouble SealevelchangeOceanAverage(Vector* Sg,SealevelMasks* masks, IssmDouble oceanarea); #endif void HydrologyEPLupdateDomainx(IssmDouble* pEplcount); Index: /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 25946) +++ /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 25947) @@ -59,5 +59,5 @@ solutioncore=&masstransport_core; break; - case SealevelriseSolutionEnum: + case SealevelchangeSolutionEnum: solutioncore=&sealevelchange_core; break; Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 25946) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 25947) @@ -54,15 +54,15 @@ void bmb_core(FemModel* femmodel); void damage_core(FemModel* femmodel); + +/*sealevel change cores:*/ void sealevelchange_core(FemModel* femmodel); void grd_core(FemModel* femmodel); +void solidearthexternal_core(FemModel* femmodel); void dynstr_core(FemModel* femmodel); -void sealevelrise_core_geometry(FemModel* femmodel); -SealevelMasks* sealevelrise_core_masks(FemModel* femmodel); -Vector* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* poceanarea); -Vector* sealevelrise_core_noneustatic(FemModel* femmodel,SealevelMasks* masks, Vector* RSLg_eustatic,IssmDouble oceanarea); -void sealevelrise_core_elastic(Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* RSLg, SealevelMasks* masks); -void sealevelrise_core_viscous(Vector** pU_gia,Vector** pN_gia,FemModel* femmodel,Vector* RSLg); -void sealevelrise_diagnostics(FemModel* femmodel,Vector* RSLg); -IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel); +void sealevelchange_geometry(FemModel* femmodel); +SealevelMasks* sealevel_masks(FemModel* femmodel); +Vector* sealevelchange_core_barystatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* poceanarea); +Vector* sealevelchange_core_sal(FemModel* femmodel,SealevelMasks* masks, Vector* RSLg_barystatic,IssmDouble oceanarea); +void sealevelchange_core_deformation(Vector** pN_radial, Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* RSLg, SealevelMasks* masks); void GetStericRate(Vector ** psteric_rate_g, FemModel* femmodel); void GetDynamicRate(Vector ** pdynamic_rate_g, FemModel* femmodel); @@ -71,4 +71,5 @@ //optimization int GradJSearch(IssmDouble* search_vector,FemModel* femmodel,int step); +IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel); //diverse Index: /issm/trunk-jpl/src/c/cores/masstransport_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 25946) +++ /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 25947) @@ -64,5 +64,5 @@ InputDuplicatex(femmodel,BaseEnum,BaseOldEnum); InputDuplicatex(femmodel,SurfaceEnum,SurfaceOldEnum); - InputDuplicatex(femmodel,SealevelriseCumDeltathicknessEnum,SealevelriseCumDeltathicknessOldEnum); + InputDuplicatex(femmodel,SealevelchangeCumDeltathicknessEnum,SealevelchangeCumDeltathicknessOldEnum); if(stabilization==4){ solutionsequence_fct(femmodel); Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 25946) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 25947) @@ -25,4 +25,6 @@ bool isgia=0; int grd=0; + int isexternal=0; + int iscoupler=0; int solution_type; @@ -32,7 +34,9 @@ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); - - /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/ - if(solution_type==SealevelriseSolutionEnum){ + femmodel->parameters->FindParam(&isexternal,SolidearthIsExternalEnum); + femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); + + /*in case we are running SealevelchangeSolutionEnum, then bypass transient settings:*/ + if(solution_type==SealevelchangeSolutionEnum){ isslr=1; isgia=1; @@ -43,5 +47,5 @@ /*Verbose: */ - if(VerboseSolution()) _printf0_(" computing sea level rise\n"); + if(VerboseSolution()) _printf0_(" computing sea level change\n"); /*Run gia core: */ @@ -55,14 +59,27 @@ /*set SLR configuration: */ - femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum); + femmodel->SetCurrentConfiguration(SealevelchangeAnalysisEnum); /*run geometry core: */ - sealevelrise_core_geometry(femmodel); + sealevelchange_geometry(femmodel); + + /*any external forcings?:*/ + if(isexternal)solidearthexternal_core(femmodel); /*Run geodetic:*/ + //if(modelid==earthid) //not sure how we proceed yet on coupling. if(grd)grd_core(femmodel); /*Run steric core for sure:*/ dynstr_core(femmodel); + + if(iscoupler){ + /*transfer sea level back to ice caps:*/ + TransferSealevel(femmodel,SealevelEnum); + TransferSealevel(femmodel,BedEnum); + + //reset cumdeltathickness to 0: + InputUpdateFromConstantx(femmodel->inputs,femmodel->elements,0.,SealevelchangeCumDeltathicknessEnum); + } /*Save results: */ @@ -71,5 +88,5 @@ char **requested_outputs = NULL; if(VerboseSolution()) _printf0_(" saving results\n"); - femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum); + femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelchangeRequestedOutputsEnum); femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} @@ -77,5 +94,5 @@ /*requested dependents: */ - if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx(); + if(solution_type==SealevelchangeSolutionEnum)femmodel->RequestedDependentsx(); /*End profiler*/ @@ -83,4 +100,89 @@ } /*}}}*/ +void solidearthexternal_core(FemModel* femmodel){ /*{{{*/ + + /*variables:*/ + Vector *bedrock = NULL; + Vector *bedrock_rate = NULL; + Vector *bedrockeast = NULL; + Vector *bedrockeast_rate = NULL; + Vector *bedrocknorth = NULL; + Vector *bedrocknorth_rate = NULL; + Vector *geoid= NULL; + Vector *geoid_rate= NULL; + int horiz=0; + int modelid=-1; + + /*parameters: */ + bool isslr=0; + int solution_type; + IssmDouble dt; + + /*Retrieve parameters:*/ + femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); + femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + + /*in case we are running SealevelchangeSolutionEnum, then bypass transient settings:*/ + if(solution_type==SealevelchangeSolutionEnum)isslr=1; + + /*Should we be here?:*/ + if(!isslr)return; + + /*Verbose: */ + if(VerboseSolution()) _printf0_(" computing external solid earth contributions\n"); + + /*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/ + GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum); + GetVectorFromInputsx(&geoid,femmodel,SealevelEnum,VertexSIdEnum); //In ISSM, Sealevel is absolute. + if(horiz){ + GetVectorFromInputsx(&bedrockeast,femmodel,BedEastEnum,VertexSIdEnum); + GetVectorFromInputsx(&bedrocknorth,femmodel,BedNorthEnum,VertexSIdEnum); + } + /*Deal with Mme: */ + if (femmodel->inputs->GetInputObjectEnum(SolidearthExternalDisplacementEastRateEnum)==DatasetInputEnum){ + /*retrieve model id: */ + femmodel->parameters->FindParam(&modelid,SolidearthExternalModelidEnum); + /*replace dataset of forcings with only one, the modelid'th:*/ + MmeToInputFromId(femmodel,modelid,SolidearthExternalDisplacementNorthRateEnum, P1Enum); + MmeToInputFromId(femmodel,modelid,SolidearthExternalDisplacementEastRateEnum, P1Enum); + MmeToInputFromId(femmodel,modelid,SolidearthExternalDisplacementUpRateEnum, P1Enum); + MmeToInputFromId(femmodel,modelid,SolidearthExternalGeoidRateEnum, P1Enum); + MmeToInputFromId(femmodel,modelid,SolidearthExternalBarystaticSeaLevelRateEnum, P1Enum); + } + + GetVectorFromInputsx(&geoid_rate,femmodel,SolidearthExternalGeoidRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&bedrock_rate,femmodel,SolidearthExternalDisplacementUpRateEnum,VertexSIdEnum); + if(horiz){ + GetVectorFromInputsx(&bedrockeast_rate,femmodel,SolidearthExternalDisplacementEastRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&bedrocknorth_rate,femmodel,SolidearthExternalDisplacementNorthRateEnum,VertexSIdEnum); + } + + /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate) * dt + steric_rate + dynamic_rate dt*/ + geoid->AXPY(geoid_rate,dt); + bedrock->AXPY(bedrock_rate,dt); + if(horiz){ + bedrockeast->AXPY(bedrockeast_rate,dt); + bedrocknorth->AXPY(bedrocknorth_rate,dt); + } + + /*update element inputs:*/ + InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,geoid,SealevelEnum,VertexSIdEnum); + if(horiz){ + InputUpdateFromVectorx(femmodel,bedrockeast,BedEastEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,bedrocknorth,BedNorthEnum,VertexSIdEnum); + } + + /*Free ressources:*/ + delete bedrock; delete bedrock_rate; + delete geoid; delete geoid_rate; + if(horiz){ + delete bedrockeast; delete bedrockeast_rate; + delete bedrocknorth; delete bedrocknorth_rate; + } +} +/*}}}*/ void grd_core(FemModel* femmodel){ /*{{{*/ @@ -89,17 +191,13 @@ /*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; + Vector *RSLg_barystatic = NULL; + Vector *U_grd = NULL; + Vector *N_grd = NULL; + Vector *U_north_grd = NULL; + Vector *U_east_grd = NULL; + Vector *bedrock = NULL; + Vector *bedrockeast = NULL; + Vector *bedrocknorth = NULL; + Vector *geoid= NULL; SealevelMasks* masks=NULL; @@ -113,7 +211,6 @@ IssmDouble dt; IssmDouble oceanarea; - int bp_compute_fingerprints=0; - - + int bp_compute_fingerprints=0; + /*Verbose: */ if(VerboseSolution()) _printf0_(" computing GRD sea level patterns\n"); @@ -122,5 +219,5 @@ femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum); - femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum); + femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum); femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); @@ -136,5 +233,5 @@ * 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 + /* in addition, if we are running solution_type SealevelchangeSolutionEnum, make sure we * run, irresepective of the time settings:*/ count=frequency; @@ -148,15 +245,15 @@ /*increment counter, or call solution core if count==frequency:*/ if (countparameters->SetParam(count,SealevelriseRunCountEnum); + count++; femmodel->parameters->SetParam(count,SealevelchangeRunCountEnum); return; } - /*call sea-level rise sub cores:*/ + /*call sea-level change sub cores:*/ if(iscoupler){ /*transfer cumulated deltathickness forcing from ice caps to earth model: */ - TransferForcing(femmodel,SealevelriseCumDeltathicknessEnum); + TransferForcing(femmodel,SealevelchangeCumDeltathicknessEnum); /*we have accumulated thicknesses, dump them in deltathcikness: */ - if(modelid==earthid)InputDuplicatex(femmodel,SealevelriseCumDeltathicknessEnum,SurfaceloadIceThicknessChangeEnum); + if(modelid==earthid)InputDuplicatex(femmodel,SealevelchangeCumDeltathicknessEnum,SurfaceloadIceThicknessChangeEnum); } @@ -165,87 +262,61 @@ /*call masks core: */ - masks=sealevelrise_core_masks(femmodel); + masks=sealevel_masks(femmodel); /*set bottom pressures:*/ SetBottomPressure(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 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) 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); - - /*transform these values into rates (as we only run this once each frequency turn:*/ - N_esa_rate=N_esa->Duplicate(); N_esa->Copy(N_esa_rate); N_esa_rate->Scale(1/(dt*frequency)); - U_esa_rate=U_esa->Duplicate(); U_esa->Copy(U_esa_rate); U_esa_rate->Scale(1/(dt*frequency)); - N_gia_rate=N_gia->Duplicate(); N_gia->Copy(N_gia_rate); N_gia_rate->Scale(1/(dt*frequency)); - U_gia_rate=U_gia->Duplicate(); U_gia->Copy(U_gia_rate); U_gia_rate->Scale(1/(dt*frequency)); - RSLg_rate=RSLg->Duplicate(); RSLg->Copy(RSLg_rate); RSLg_rate->Scale(1/(dt*frequency)); - - /*get some results into elements:{{{*/ - InputUpdateFromVectorx(femmodel,U_esa_rate,SealevelUEsaRateEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,N_esa_rate,SealevelNEsaRateEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,U_gia_rate,UGiaRateEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,N_gia_rate,NGiaRateEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,RSLg_rate,SealevelRSLRateEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,U_esa,SealevelUEsaEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,N_esa,SealevelNEsaEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,U_gia,UGiaEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,N_gia,NGiaEnum,VertexSIdEnum); + /*call barystatic core (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */ + RSLg_barystatic=sealevelchange_core_barystatic(femmodel,masks,&oceanarea); + + /*call self attraction and loading module (ocean loading tems - 2nd and 5th terms on the RHS of Farrel and Clark) */ + RSLg=sealevelchange_core_sal(femmodel,masks,RSLg_barystatic,oceanarea); + + /*compute bedrock motion and derive geoid: */ + sealevelchange_core_deformation(&N_grd,&U_grd,&U_north_grd,&U_east_grd,femmodel,RSLg,masks); + + /*Update bedrock motion and geoid:*/ + GetVectorFromInputsx(&geoid,femmodel,SealevelEnum,VertexSIdEnum); + GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum); + if(horiz){ + GetVectorFromInputsx(&bedrockeast,femmodel,BedEastEnum,VertexSIdEnum); + GetVectorFromInputsx(&bedrocknorth,femmodel,BedNorthEnum,VertexSIdEnum); + } + + geoid->AXPY(N_grd,1); + bedrock->AXPY(U_grd,1); + if(horiz){ + bedrockeast->AXPY(U_east_grd,1); + bedrocknorth->AXPY(U_north_grd,1); + } + + /*get some of the updates into elements:*/ + InputUpdateFromVectorx(femmodel,U_grd,SealevelUEsaEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,N_grd,SealevelNEsaEnum,VertexSIdEnum); InputUpdateFromVectorx(femmodel,RSLg,SealevelRSLEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,RSLg_eustatic,SealevelRSLEustaticEnum,VertexSIdEnum); - + InputUpdateFromVectorx(femmodel,RSLg_barystatic,SealevelRSLBarystaticEnum,VertexSIdEnum); if (horiz){ - InputUpdateFromVectorx(femmodel,U_north_esa,SealevelUNorthEsaEnum,VertexSIdEnum); // north motion - InputUpdateFromVectorx(femmodel,U_east_esa,SealevelUEastEsaEnum,VertexSIdEnum); // east motion - } /*}}}*/ - } - - if(iscoupler){ - /*transfer sea level back to ice caps:*/ - TransferSealevel(femmodel,SealevelNEsaRateEnum); - TransferSealevel(femmodel,NGiaRateEnum); - TransferSealevel(femmodel,SealevelUEsaRateEnum); - TransferSealevel(femmodel,UGiaRateEnum); - - //reset cumdeltathickness to 0: - InputUpdateFromConstantx(femmodel->inputs,femmodel->elements,0.,SealevelriseCumDeltathicknessEnum); + InputUpdateFromVectorx(femmodel,U_north_grd,SealevelUNorthEsaEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,U_east_grd,SealevelUEastEsaEnum,VertexSIdEnum); + } } /*reset counter to 1:*/ - femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter. + femmodel->parameters->SetParam(1,SealevelchangeRunCountEnum); //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; - + delete RSLg_barystatic; + delete U_grd; + delete N_grd; + delete bedrock; + delete geoid; 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; + delete U_north_grd; + delete U_east_grd; + delete bedrockeast; + delete bedrocknorth; + } + delete masks; /*}}}*/ @@ -255,12 +326,7 @@ /*variables:*/ - Vector *bedrock = NULL; - Vector *SL = NULL; + Vector *sealevel = NULL; Vector *steric_rate_g = NULL; Vector *dynamic_rate_g = NULL; - Vector *U_esa_rate= NULL; - Vector *N_esa_rate= NULL; - Vector *U_gia_rate= NULL; - Vector *N_gia_rate= NULL; /*parameters: */ @@ -268,7 +334,7 @@ int solution_type; IssmDouble dt; - int grd=0; int step; IssmDouble time; + IssmDouble cumgmtslr=0; IssmDouble cumbslr=0; @@ -277,76 +343,57 @@ /*Retrieve parameters:*/ + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); + + /*in case we are running SealevelchangeSolutionEnum, then bypass transient settings:*/ + if(solution_type==SealevelchangeSolutionEnum)isslr=1; + + /*Should we be here?:*/ + if(!isslr)return; + + /*Verbose: */ + if(VerboseSolution()) _printf0_(" computing steric sea level change\n"); + + /*Retrieve sealevel and add steric + dynamic rates:*/ + GetVectorFromInputsx(&sealevel,femmodel,SealevelEnum,VertexSIdEnum); + GetStericRate(&steric_rate_g,femmodel); + GetDynamicRate(&dynamic_rate_g,femmodel); + + /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate) * dt + steric_rate + dynamic_rate dt*/ femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); - femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); + sealevel->AXPY(steric_rate_g,dt); + sealevel->AXPY(dynamic_rate_g,dt); + + /*cumulate thermal steric rate:*/ femmodel->parameters->FindParam(&cumgmtslr,CumGmtslrEnum); femmodel->parameters->FindParam(&cumbslr,CumBslrEnum); - femmodel->parameters->FindParam(&step,StepEnum); - femmodel->parameters->FindParam(&time,TimeEnum); - - /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/ - if(solution_type==SealevelriseSolutionEnum)isslr=1; - - /*Should we be here?:*/ - if(!isslr)return; - - /*Verbose: */ - if(VerboseSolution()) _printf0_(" computing steric sea level rise\n"); - - /*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/ - GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum); - GetVectorFromInputsx(&SL,femmodel,SealevelEnum,VertexSIdEnum); - GetStericRate(&steric_rate_g,femmodel); - GetDynamicRate(&dynamic_rate_g,femmodel); - if(grd){ - GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum); - GetVectorFromInputsx(&U_gia_rate,femmodel,UGiaRateEnum,VertexSIdEnum); - GetVectorFromInputsx(&N_esa_rate,femmodel,SealevelNEsaRateEnum,VertexSIdEnum); - GetVectorFromInputsx(&N_gia_rate,femmodel,NGiaRateEnum,VertexSIdEnum); - } - - /*Save steric rate, and cumulate it through time, cumulate with barystatic sea level rise, and cumulate through time:*/ + gmtslr=steric_rate_g->Max()*dt; cumgmtslr+=gmtslr; cumgmslr=cumbslr+cumgmtslr; + + femmodel->parameters->SetParam(cumgmtslr,CumGmtslrEnum); + femmodel->parameters->SetParam(cumgmslr,CumGmslrEnum); + + /*Outputs some metrics:*/ + femmodel->parameters->FindParam(&step,StepEnum); + femmodel->parameters->FindParam(&time,TimeEnum); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,GmtslrEnum,gmtslr,step,time)); femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,CumGmtslrEnum,cumgmtslr,step,time)); femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,CumGmslrEnum,cumgmslr,step,time)); - femmodel->parameters->SetParam(cumgmtslr,CumGmtslrEnum); - - /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate) * dt + steric_rate + dynamic_rate dt*/ - if(grd){ - SL->AXPY(N_gia_rate,dt); - SL->AXPY(N_esa_rate,dt); - } - SL->AXPY(steric_rate_g,dt); - SL->AXPY(dynamic_rate_g,dt); - - /*compute new bedrock position: */ - if(grd){ - 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); + InputUpdateFromVectorx(femmodel,sealevel,SealevelEnum,VertexSIdEnum); /*Free ressources:*/ - delete bedrock; - delete SL; + delete sealevel; delete steric_rate_g; delete dynamic_rate_g; - if(grd){ - delete U_esa_rate; - delete U_gia_rate; - delete N_esa_rate; - delete N_gia_rate; - } } /*}}}*/ -SealevelMasks* sealevelrise_core_masks(FemModel* femmodel) { /*{{{*/ +//Geometry: +SealevelMasks* sealevel_masks(FemModel* femmodel) { /*{{{*/ if(VerboseSolution()) _printf0_(" computing sea level masks\n"); @@ -363,5 +410,5 @@ return masks; }/*}}}*/ -void sealevelrise_core_geometry(FemModel* femmodel) { /*{{{*/ +void sealevelchange_geometry(FemModel* femmodel) { /*{{{*/ /*Geometry core where we compute indices into tables pre computed in the SealevelRiseAnalysis: */ @@ -381,5 +428,5 @@ /*retrieve parameters:*/ femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); - femmodel->parameters->FindParam(&geometrydone,SealevelriseGeometryDoneEnum); + femmodel->parameters->FindParam(&geometrydone,SealevelchangeGeometryDoneEnum); if(geometrydone){ @@ -396,8 +443,8 @@ - /*Run sealevelrise geometry routine in elements:*/ + /*Run sealevel geometry routine in elements:*/ for(Object* & object : femmodel->elements->objects){ Element* element=xDynamicCast(object); - element->SealevelriseGeometry(latitude,longitude,radius,xx,yy,zz); + element->SealevelchangeGeometry(latitude,longitude,radius,xx,yy,zz); } @@ -413,11 +460,13 @@ /*Record the fact that we ran this module already: */ - femmodel->parameters->SetParam(true,SealevelriseGeometryDoneEnum); + femmodel->parameters->SetParam(true,SealevelchangeGeometryDoneEnum); }/*}}}*/ -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)*/ + +//GRD: +Vector* sealevelchange_core_barystatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* poceanarea){ /*{{{*/ + + /*Barystatic core of the SLR solution (terms that are constant with respect to sea-level)*/ Vector *RSLgi = NULL; @@ -464,8 +513,8 @@ /*call the bslr main module: */ - femmodel->SealevelriseEustatic(RSLgi,&oceanarea,&bslr, &bslrice, &bslrhydro, &bslrice_partition, &bslrhydro_partition,masks); //this computes + femmodel->SealevelchangeBarystatic(RSLgi,&oceanarea,&bslr, &bslrice, &bslrhydro, &bslrice_partition, &bslrhydro_partition,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); + RSLgi_oceanaverage=femmodel->SealevelchangeOceanAverage(RSLgi,masks, oceanarea); /*RSLg is the sum of the pure bslr component (term 3) and the contribution from the perturbation to the graviation potential due to the @@ -509,8 +558,8 @@ 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* sealevelchange_core_sal(FemModel* femmodel, SealevelMasks* masks, Vector* RSLg_barystatic,IssmDouble oceanarea){ /*{{{*/ + + /*this core computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5. + sal core of the SLR solution */ Vector *RSLg = NULL; @@ -549,5 +598,5 @@ RSLg = new Vector(gsize); RSLg->Assemble(); - RSLg_eustatic->Copy(RSLg); //first initialize RSLg with the eustatic component computed in sealevelrise_core_eustatic. + RSLg_barystatic->Copy(RSLg); //first initialize RSLg with the barystatic component computed in sealevelchange_core_barystatic. RSLg_old = new Vector(gsize); @@ -560,5 +609,5 @@ for(;;){ - //save pointer to old sea level rise + //save pointer to old sea level delete RSLg_old; RSLg_old=RSLg; @@ -567,6 +616,6 @@ RSLgo = new Vector(gsize); RSLgo->Assemble(); - /*call the non eustatic module: */ - femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, masks, verboseconvolution); + /*call the sal module: */ + femmodel->SealevelchangeSal(RSLgo, RSLg_old, masks, verboseconvolution); /*assemble solution vector: */ @@ -577,5 +626,5 @@ /*call rotational feedback module: */ RSLgo_rot = new Vector(gsize); RSLgo_rot->Assemble(); - femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, masks); + femmodel->SealevelchangeRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, masks); RSLgo_rot->Assemble(); @@ -589,8 +638,8 @@ /*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); + RSLgo_oceanaverage=femmodel->SealevelchangeOceanAverage(RSLgo,masks, oceanarea); + + /*RSLg is the sum of the barystatic term, and the ocean terms: */ + RSLg_barystatic->Copy(RSLg); RSLg->AXPY(RSLgo,1); RSLg->Shift(-RSLgo_oceanaverage); @@ -625,9 +674,10 @@ 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; +void sealevelchange_core_deformation(Vector** pN_grd, Vector** pU_grd, Vector** pU_north_grd,Vector** pU_east_grd,FemModel* femmodel,Vector* RSLg, SealevelMasks* masks){ /*{{{*/ + + Vector *U_grd = NULL; + Vector *N_grd = NULL; + Vector *U_north_grd = NULL; + Vector *U_east_grd = NULL; /*parameters: */ @@ -652,10 +702,11 @@ /*intialize vectors:*/ - U_esa = new Vector(gsize); + U_grd = new Vector(gsize); + N_grd = new Vector(gsize); if (horiz){ - U_north_esa = new Vector(gsize); - U_east_esa = new Vector(gsize); - } - + U_north_grd = new Vector(gsize); + U_east_grd = new Vector(gsize); + } + /*retrieve geometric information: */ VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); @@ -663,11 +714,12 @@ /*call the elastic main modlule:*/ - femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg, masks); + femmodel->SealevelchangeDeformation(N_grd, U_grd,U_north_grd,U_east_grd,RSLg, masks); /*Assign output pointers:*/ - *pU_esa=U_esa; + *pU_grd=U_grd; + *pN_grd=N_grd; if(horiz){ - *pU_east_esa=U_east_esa; - *pU_north_esa=U_north_esa; + *pU_east_grd=U_east_grd; + *pU_north_grd=U_north_grd; } @@ -681,40 +733,6 @@ } /*}}}*/ -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,SolidearthSettingsRunFrequencyEnum); - femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - - /*recover GIA rates:*/ - GetVectorFromInputsx(&U_gia,femmodel,UGiaRateEnum,VertexSIdEnum); - GetVectorFromInputsx(&N_gia,femmodel,NGiaRateEnum,VertexSIdEnum); - - /*we just loaded rates, that's not what's being asked, scale by time:*/ - U_gia->Scale(frequency*dt); - N_gia->Scale(frequency*dt); - - /*Assign output pointers:*/ - *pU_gia=U_gia; - *pN_gia=N_gia; - -} -/*}}}*/ -void sealevelrise_diagnostics(FemModel* femmodel,Vector* RSLg){ /*{{{*/ - - /*compute spherical harmonics deg 1 and deg 2 coefficeints:*/ - -} -/*}}}*/ + +//Ocean: void GetDynamicRate(Vector ** pdynamic_rate_g, FemModel* femmodel){ /*{{{*/ @@ -837,5 +855,5 @@ /*}}}*/ -/*support routines:*/ +/*Support routines:*/ void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/ @@ -913,5 +931,5 @@ /*Retrieve transition vectors, used to plug from each ice cap into the global forcing:*/ - femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum); + femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelchangeTransitionsEnum); if(ntransitions!=earthid)_error_("TransferForcing error message: number of transition vectors is not equal to the number of icecaps!"); @@ -1019,5 +1037,5 @@ /*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); + femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelchangeTransitionsEnum); if(ntransitions!=earthid)_error_("TransferSeaLevel error message: number of transition vectors is not equal to the number of icecaps!"); @@ -1087,5 +1105,5 @@ /*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.*/ + * the thickness variations from SealevelchangeSpcthicknessEnum.*/ /*No mass transport module was called, so we are just going to retrieve the geometry thickness @@ -1095,6 +1113,6 @@ nv=femmodel->vertices->NumberOfVertices(); - GetVectorFromInputsx(&newthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum); - GetVectorFromInputsx(&oldthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum,time-dt);*/ + GetVectorFromInputsx(&newthickness,femmodel,SealevelchangeSpcthicknessEnum,VertexSIdEnum); + GetVectorFromInputsx(&oldthickness,femmodel,SealevelchangeSpcthicknessEnum,VertexSIdEnum,time-dt);*/ /*compute deltathickness: */ @@ -1106,7 +1124,7 @@ /*add to cumulated delta thickness: */ - /*GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); + /*GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelchangeCumDeltathicknessEnum,VertexSIdEnum); cumdeltathickness->AXPY(deltathickness,1); - InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);*/ + InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelchangeCumDeltathicknessEnum,VertexSIdEnum);*/ /*free ressources:*/ Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 25946) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 25947) @@ -62,5 +62,5 @@ #endif - if(isslr) sealevelrise_core_geometry(femmodel); + if(isslr) sealevelchange_geometry(femmodel); while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime. @@ -235,7 +235,7 @@ if(isesa) esa_core(femmodel); - /*Sea level rise: */ + /*Sea level change: */ if(isslr){ - if(VerboseSolution()) _printf0_(" computing sea level rise\n"); + if(VerboseSolution()) _printf0_(" computing sea level change\n"); sealevelchange_core(femmodel); } @@ -275,5 +275,5 @@ femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); femmodel->parameters->FindParam(&checkpoint_frequency,SettingsCheckpointFrequencyEnum); _assert_(checkpoint_frequency>0); - if(isslr) sealevelrise_core_geometry(femmodel); + if(isslr) sealevelchange_geometry(femmodel); std::vector time_all; Index: /issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp (revision 25946) +++ /issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp (revision 25947) @@ -160,5 +160,5 @@ /*Here, we put all the code that cannot be handled any other place: */ - if (strncmp(root,"SurfaceloadModelid",18)==0){ + if (strncmp(root,"SurfaceloadModelid",18)==0){ //surface load in solid earth class {{{ if(VerboseQmu()){ @@ -178,5 +178,29 @@ MmeToInput(femmodel,distributed_values,variable_partition,npart,MaskOceanLevelsetEnum, P1Enum); - } + } /*}}}*/ + if (strncmp(root,"SolidearthExternalModelid",18)==0){ //external solid earth solution in solid earth class {{{ + + if(VerboseQmu()){ + _printf0_(" SolidearthExternalModelid MME, with ids: "); + for (int i=0;iinputs->GetInputObjectEnum(SolidearthExternalDisplacementEastRateEnum)==DatasetInputEnum) + MmeToInput(femmodel,distributed_values,variable_partition,npart,SolidearthExternalDisplacementEastRateEnum, P1Enum); + + if (femmodel->inputs->GetInputObjectEnum(SolidearthExternalDisplacementUpRateEnum)==DatasetInputEnum) + MmeToInput(femmodel,distributed_values,variable_partition,npart,SolidearthExternalDisplacementUpRateEnum, P1Enum); + + if (femmodel->inputs->GetInputObjectEnum(SolidearthExternalDisplacementNorthRateEnum)==DatasetInputEnum) + MmeToInput(femmodel,distributed_values,variable_partition,npart,SolidearthExternalDisplacementNorthRateEnum, P1Enum); + + if (femmodel->inputs->GetInputObjectEnum(SolidearthExternalGeoidRateEnum)==DatasetInputEnum) + MmeToInput(femmodel,distributed_values,variable_partition,npart,SolidearthExternalGeoidRateEnum, P1Enum); + + if (femmodel->inputs->GetInputObjectEnum(SolidearthExternalBarystaticSeaLevelRateEnum)==DatasetInputEnum) + MmeToInput(femmodel,distributed_values,variable_partition,npart,SolidearthExternalBarystaticSeaLevelRateEnum, P1Enum); + } /*}}}*/ + else _error_("InputUpdateSpecialtyCode error message: " << root << " not supported yet!"); @@ -241,4 +265,59 @@ transientinput2->Configure(femmodel->parameters); } //}}} +void MmeToInputFromId(FemModel* femmodel,int id, int rootenum, int interpolationenum){ /*{{{*/ + + TransientInput* transientinput = NULL; + TransientInput* transientinput2 = NULL; + Tria* element = NULL; + IssmDouble value; + IssmDouble* values = NULL; + IssmDouble* times = NULL; + int N; + + /*find thickness dataset: */ + DatasetInput* datasetinput = femmodel->inputs->GetDatasetInput(rootenum); + + /*Initialize new transient input: */ + transientinput = datasetinput->GetTransientInputByOffset(0); _assert_(transientinput); + transientinput->GetAllTimes(×,&N); + femmodel->inputs->SetTransientInput(DummyEnum,times,N); + transientinput2 = femmodel->inputs->GetTransientInput(DummyEnum); + + for(Object* & object : femmodel->elements->objects){ + Tria* element=xDynamicCast(object); + + /*recover the right field from the mme: */ + transientinput = datasetinput->GetTransientInputByOffset(id); _assert_(transientinput); + + /*copy values from the transientinput to the final transientinput2: */ + for (int j=0;jGetTriaInput(j); + element->InputServe(tria_input); + if(interpolationenum==P0Enum){ + value=tria_input->element_values[0]; + transientinput2->AddTriaTimeInput( j,1,&(element->lid),&value,P0Enum); + } + else if(interpolationenum==P1Enum){ + + /*Get values and lid list*/ + const int numvertices = element->GetNumberOfVertices(); + int *vertexlids = xNew(numvertices); + int *vertexsids = xNew(numvertices); + + /*Recover vertices ids needed to initialize inputs*/ + element->GetVerticesLidList(&vertexlids[0]); + element->GetVerticesSidList(&vertexsids[0]); + values=tria_input->element_values; + transientinput2->AddTriaTimeInput( j,numvertices,vertexlids,values,P1Enum); + } + } + } + + /*wipe out existing SurfaceloadIceThicknessChangeEnum dataset:*/ + femmodel->inputs->ChangeEnum(DummyEnum,rootenum); + + //reconfigure: + transientinput2->Configure(femmodel->parameters); +} //}}} void InputScaleFromDakotax(FemModel* femmodel,IssmDouble* distributed_values,IssmDouble* partition, int npart, int nt, int name){ /*{{{*/ Index: /issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.h =================================================================== --- /issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.h (revision 25946) +++ /issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.h (revision 25947) @@ -11,4 +11,5 @@ void InputUpdateSpecialtyCode(FemModel* femmodel,IssmDouble* distributed_values,IssmDouble* variable_partition,int npart,char* root); void MmeToInput(FemModel* femmodel,IssmDouble* distributed_values,IssmDouble* variable_partition,int npart,int rootenum, int interpolationenum); +void MmeToInputFromId(FemModel* femmodel,int id, int rootenum, int interpolationenum); void InputScaleFromDakotax(FemModel* femmodel,IssmDouble* distributed_values,IssmDouble* partition, int npart, int nt, int name); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 25946) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 25947) @@ -63,5 +63,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.parameterization",FrontalForcingsParamEnum)); - parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1)); + parameters->AddObject(new IntParam(SealevelchangeRunCountEnum,1)); {/*This is specific to ice...*/ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25946) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25947) @@ -138,4 +138,8 @@ syn keyword cConstant DslModelidEnum syn keyword cConstant DslNummodelsEnum +syn keyword cConstant SolidearthIsExternalEnum +syn keyword cConstant SolidearthExternalNatureEnum +syn keyword cConstant SolidearthExternalModelidEnum +syn keyword cConstant SolidearthExternalNummodelsEnum syn keyword cConstant DslComputeFingerprintsEnum syn keyword cConstant EarthIdEnum @@ -338,5 +342,5 @@ syn keyword cConstant RotationalAngularVelocityEnum syn keyword cConstant SolidearthSettingsElasticEnum -syn keyword cConstant SealevelriseGeometryDoneEnum +syn keyword cConstant SealevelchangeGeometryDoneEnum syn keyword cConstant RotationalEquatorialMoiEnum syn keyword cConstant TidalLoveHEnum @@ -347,11 +351,11 @@ syn keyword cConstant LoadLoveKEnum syn keyword cConstant LoadLoveLEnum -syn keyword cConstant SealevelriseGRigidEnum -syn keyword cConstant SealevelriseGElasticEnum +syn keyword cConstant SealevelchangeGRigidEnum +syn keyword cConstant SealevelchangeGElasticEnum syn keyword cConstant SolidearthSettingsComputesealevelchangeEnum syn keyword cConstant SolidearthSettingsGRDEnum syn keyword cConstant SolidearthSettingsGlfractionEnum syn keyword cConstant SolidearthSettingsRunFrequencyEnum -syn keyword cConstant SealevelriseHElasticEnum +syn keyword cConstant SealevelchangeHElasticEnum syn keyword cConstant SolidearthSettingsHorizEnum syn keyword cConstant SolidearthSettingsMaxiterEnum @@ -359,10 +363,10 @@ syn keyword cConstant RotationalPolarMoiEnum syn keyword cConstant SolidearthSettingsReltolEnum -syn keyword cConstant SealevelriseRequestedOutputsEnum +syn keyword cConstant SealevelchangeRequestedOutputsEnum syn keyword cConstant SolidearthSettingsRigidEnum syn keyword cConstant SolidearthSettingsRotationEnum -syn keyword cConstant SealevelriseRunCountEnum -syn keyword cConstant SealevelriseTransitionsEnum -syn keyword cConstant SealevelriseUElasticEnum +syn keyword cConstant SealevelchangeRunCountEnum +syn keyword cConstant SealevelchangeTransitionsEnum +syn keyword cConstant SealevelchangeUElasticEnum syn keyword cConstant SettingsIoGatherEnum syn keyword cConstant SettingsNumResultsOnNodesEnum @@ -532,4 +536,6 @@ syn keyword cConstant BaseSlopeYEnum syn keyword cConstant BedEnum +syn keyword cConstant BedEastEnum +syn keyword cConstant BedNorthEnum syn keyword cConstant BedSlopeXEnum syn keyword cConstant BedSlopeYEnum @@ -710,6 +716,6 @@ syn keyword cConstant RheologyBbarAbsGradientEnum syn keyword cConstant SealevelEnum -syn keyword cConstant SealevelEustaticMaskEnum -syn keyword cConstant SealevelEustaticOceanMaskEnum +syn keyword cConstant SealevelBarystaticMaskEnum +syn keyword cConstant SealevelBarystaticOceanMaskEnum syn keyword cConstant SealevelNEsaEnum syn keyword cConstant SealevelNEsaRateEnum @@ -720,6 +726,5 @@ syn keyword cConstant BslrRateEnum syn keyword cConstant GmtslrEnum -syn keyword cConstant SealevelRSLEustaticEnum -syn keyword cConstant SealevelRSLEustaticRateEnum +syn keyword cConstant SealevelRSLBarystaticEnum syn keyword cConstant SealevelRSLRateEnum syn keyword cConstant SealevelUEastEsaEnum @@ -727,14 +732,14 @@ syn keyword cConstant SealevelUEsaRateEnum syn keyword cConstant SealevelUNorthEsaEnum -syn keyword cConstant SealevelriseCumDeltathicknessEnum -syn keyword cConstant SealevelriseCumDeltathicknessOldEnum +syn keyword cConstant SealevelchangeCumDeltathicknessEnum +syn keyword cConstant SealevelchangeCumDeltathicknessOldEnum syn keyword cConstant SurfaceloadOtherEnum syn keyword cConstant SurfaceloadIceThicknessChangeEnum syn keyword cConstant SurfaceloadWaterHeightChangeEnum -syn keyword cConstant SealevelriseIndicesEnum -syn keyword cConstant SealevelriseGEnum -syn keyword cConstant SealevelriseGUEnum -syn keyword cConstant SealevelriseGEEnum -syn keyword cConstant SealevelriseGNEnum +syn keyword cConstant SealevelchangeIndicesEnum +syn keyword cConstant SealevelchangeGEnum +syn keyword cConstant SealevelchangeGUEnum +syn keyword cConstant SealevelchangeGEEnum +syn keyword cConstant SealevelchangeGNEnum syn keyword cConstant SedimentHeadEnum syn keyword cConstant SedimentHeadOldEnum @@ -844,4 +849,9 @@ syn keyword cConstant SmbZTopEnum syn keyword cConstant SmbZYEnum +syn keyword cConstant SolidearthExternalDisplacementEastRateEnum +syn keyword cConstant SolidearthExternalDisplacementNorthRateEnum +syn keyword cConstant SolidearthExternalDisplacementUpRateEnum +syn keyword cConstant SolidearthExternalGeoidRateEnum +syn keyword cConstant SolidearthExternalBarystaticSeaLevelRateEnum syn keyword cConstant StrainRateeffectiveEnum syn keyword cConstant StrainRateparallelEnum @@ -1303,6 +1313,6 @@ syn keyword cConstant SealevelNmotionEnum syn keyword cConstant SealevelUmotionEnum -syn keyword cConstant SealevelriseAnalysisEnum -syn keyword cConstant SealevelriseSolutionEnum +syn keyword cConstant SealevelchangeAnalysisEnum +syn keyword cConstant SealevelchangeSolutionEnum syn keyword cConstant SegEnum syn keyword cConstant SegInputEnum @@ -1399,5 +1409,4 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel -syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1406,6 +1415,6 @@ syn keyword cType ControlInput syn keyword cType Covertree +syn keyword cType DataSetParam syn keyword cType DatasetInput -syn keyword cType DataSetParam syn keyword cType Definition syn keyword cType DependentObject @@ -1420,6 +1429,6 @@ syn keyword cType ElementInput syn keyword cType ElementMatrix +syn keyword cType ElementVector syn keyword cType Elements -syn keyword cType ElementVector syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1428,10 +1437,9 @@ syn keyword cType Friction syn keyword cType Gauss -syn keyword cType GaussianVariogram -syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria +syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1448,5 +1456,4 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface -syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1459,5 +1466,4 @@ syn keyword cType Matice syn keyword cType Matlitho -syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1472,6 +1478,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType OptionUtilities syn keyword cType Options -syn keyword cType OptionUtilities syn keyword cType Param syn keyword cType Parameters @@ -1487,11 +1493,11 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType RiftStruct syn keyword cType Riftfront -syn keyword cType RiftStruct syn keyword cType SealevelMasks syn keyword cType Seg syn keyword cType SegInput +syn keyword cType SegRef syn keyword cType Segment -syn keyword cType SegRef syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1512,4 +1518,8 @@ syn keyword cType Vertex syn keyword cType Vertices +syn keyword cType classes +syn keyword cType gaussobjects +syn keyword cType krigingobjects +syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis @@ -1530,6 +1540,6 @@ syn keyword cType FreeSurfaceBaseAnalysis syn keyword cType FreeSurfaceTopAnalysis +syn keyword cType GLheightadvectionAnalysis syn keyword cType GiaAnalysis -syn keyword cType GLheightadvectionAnalysis syn keyword cType HydrologyDCEfficientAnalysis syn keyword cType HydrologyDCInefficientAnalysis @@ -1544,5 +1554,5 @@ syn keyword cType MasstransportAnalysis syn keyword cType MeltingAnalysis -syn keyword cType SealevelriseAnalysis +syn keyword cType SealevelchangeAnalysis syn keyword cType SmbAnalysis syn keyword cType SmoothAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25946) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25947) @@ -132,4 +132,8 @@ DslModelidEnum, DslNummodelsEnum, + SolidearthIsExternalEnum, + SolidearthExternalNatureEnum, + SolidearthExternalModelidEnum, + SolidearthExternalNummodelsEnum, DslComputeFingerprintsEnum, EarthIdEnum, @@ -332,5 +336,5 @@ RotationalAngularVelocityEnum, SolidearthSettingsElasticEnum, - SealevelriseGeometryDoneEnum, + SealevelchangeGeometryDoneEnum, RotationalEquatorialMoiEnum, TidalLoveHEnum, @@ -341,11 +345,11 @@ LoadLoveKEnum, LoadLoveLEnum, - SealevelriseGRigidEnum, - SealevelriseGElasticEnum, + SealevelchangeGRigidEnum, + SealevelchangeGElasticEnum, SolidearthSettingsComputesealevelchangeEnum, SolidearthSettingsGRDEnum, SolidearthSettingsGlfractionEnum, SolidearthSettingsRunFrequencyEnum, - SealevelriseHElasticEnum, + SealevelchangeHElasticEnum, SolidearthSettingsHorizEnum, SolidearthSettingsMaxiterEnum, @@ -353,10 +357,10 @@ RotationalPolarMoiEnum, SolidearthSettingsReltolEnum, - SealevelriseRequestedOutputsEnum, + SealevelchangeRequestedOutputsEnum, SolidearthSettingsRigidEnum, SolidearthSettingsRotationEnum, - SealevelriseRunCountEnum, - SealevelriseTransitionsEnum, - SealevelriseUElasticEnum, + SealevelchangeRunCountEnum, + SealevelchangeTransitionsEnum, + SealevelchangeUElasticEnum, SettingsIoGatherEnum, SettingsNumResultsOnNodesEnum, @@ -528,4 +532,6 @@ BaseSlopeYEnum, BedEnum, + BedEastEnum, + BedNorthEnum, BedSlopeXEnum, BedSlopeYEnum, @@ -706,6 +712,6 @@ RheologyBbarAbsGradientEnum, SealevelEnum, - SealevelEustaticMaskEnum, - SealevelEustaticOceanMaskEnum, + SealevelBarystaticMaskEnum, + SealevelBarystaticOceanMaskEnum, SealevelNEsaEnum, SealevelNEsaRateEnum, @@ -716,6 +722,5 @@ BslrRateEnum, GmtslrEnum, - SealevelRSLEustaticEnum, - SealevelRSLEustaticRateEnum, + SealevelRSLBarystaticEnum, SealevelRSLRateEnum, SealevelUEastEsaEnum, @@ -723,14 +728,14 @@ SealevelUEsaRateEnum, SealevelUNorthEsaEnum, - SealevelriseCumDeltathicknessEnum, - SealevelriseCumDeltathicknessOldEnum, + SealevelchangeCumDeltathicknessEnum, + SealevelchangeCumDeltathicknessOldEnum, SurfaceloadOtherEnum, SurfaceloadIceThicknessChangeEnum, SurfaceloadWaterHeightChangeEnum, - SealevelriseIndicesEnum, - SealevelriseGEnum, - SealevelriseGUEnum, - SealevelriseGEEnum, - SealevelriseGNEnum, + SealevelchangeIndicesEnum, + SealevelchangeGEnum, + SealevelchangeGUEnum, + SealevelchangeGEEnum, + SealevelchangeGNEnum, SedimentHeadEnum, SedimentHeadOldEnum, @@ -841,4 +846,9 @@ SmbZTopEnum, SmbZYEnum, + SolidearthExternalDisplacementEastRateEnum, + SolidearthExternalDisplacementNorthRateEnum, + SolidearthExternalDisplacementUpRateEnum, + SolidearthExternalGeoidRateEnum, + SolidearthExternalBarystaticSeaLevelRateEnum, StrainRateeffectiveEnum, StrainRateparallelEnum, @@ -1302,6 +1312,6 @@ SealevelNmotionEnum, SealevelUmotionEnum, - SealevelriseAnalysisEnum, - SealevelriseSolutionEnum, + SealevelchangeAnalysisEnum, + SealevelchangeSolutionEnum, SegEnum, SegInputEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25946) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25947) @@ -140,4 +140,8 @@ case DslModelidEnum : return "DslModelid"; case DslNummodelsEnum : return "DslNummodels"; + case SolidearthIsExternalEnum : return "SolidearthIsExternal"; + case SolidearthExternalNatureEnum : return "SolidearthExternalNature"; + case SolidearthExternalModelidEnum : return "SolidearthExternalModelid"; + case SolidearthExternalNummodelsEnum : return "SolidearthExternalNummodels"; case DslComputeFingerprintsEnum : return "DslComputeFingerprints"; case EarthIdEnum : return "EarthId"; @@ -340,5 +344,5 @@ case RotationalAngularVelocityEnum : return "RotationalAngularVelocity"; case SolidearthSettingsElasticEnum : return "SolidearthSettingsElastic"; - case SealevelriseGeometryDoneEnum : return "SealevelriseGeometryDone"; + case SealevelchangeGeometryDoneEnum : return "SealevelchangeGeometryDone"; case RotationalEquatorialMoiEnum : return "RotationalEquatorialMoi"; case TidalLoveHEnum : return "TidalLoveH"; @@ -349,11 +353,11 @@ case LoadLoveKEnum : return "LoadLoveK"; case LoadLoveLEnum : return "LoadLoveL"; - case SealevelriseGRigidEnum : return "SealevelriseGRigid"; - case SealevelriseGElasticEnum : return "SealevelriseGElastic"; + case SealevelchangeGRigidEnum : return "SealevelchangeGRigid"; + case SealevelchangeGElasticEnum : return "SealevelchangeGElastic"; case SolidearthSettingsComputesealevelchangeEnum : return "SolidearthSettingsComputesealevelchange"; case SolidearthSettingsGRDEnum : return "SolidearthSettingsGRD"; case SolidearthSettingsGlfractionEnum : return "SolidearthSettingsGlfraction"; case SolidearthSettingsRunFrequencyEnum : return "SolidearthSettingsRunFrequency"; - case SealevelriseHElasticEnum : return "SealevelriseHElastic"; + case SealevelchangeHElasticEnum : return "SealevelchangeHElastic"; case SolidearthSettingsHorizEnum : return "SolidearthSettingsHoriz"; case SolidearthSettingsMaxiterEnum : return "SolidearthSettingsMaxiter"; @@ -361,10 +365,10 @@ case RotationalPolarMoiEnum : return "RotationalPolarMoi"; case SolidearthSettingsReltolEnum : return "SolidearthSettingsReltol"; - case SealevelriseRequestedOutputsEnum : return "SealevelriseRequestedOutputs"; + case SealevelchangeRequestedOutputsEnum : return "SealevelchangeRequestedOutputs"; case SolidearthSettingsRigidEnum : return "SolidearthSettingsRigid"; case SolidearthSettingsRotationEnum : return "SolidearthSettingsRotation"; - case SealevelriseRunCountEnum : return "SealevelriseRunCount"; - case SealevelriseTransitionsEnum : return "SealevelriseTransitions"; - case SealevelriseUElasticEnum : return "SealevelriseUElastic"; + case SealevelchangeRunCountEnum : return "SealevelchangeRunCount"; + case SealevelchangeTransitionsEnum : return "SealevelchangeTransitions"; + case SealevelchangeUElasticEnum : return "SealevelchangeUElastic"; case SettingsIoGatherEnum : return "SettingsIoGather"; case SettingsNumResultsOnNodesEnum : return "SettingsNumResultsOnNodes"; @@ -534,4 +538,6 @@ case BaseSlopeYEnum : return "BaseSlopeY"; case BedEnum : return "Bed"; + case BedEastEnum : return "BedEast"; + case BedNorthEnum : return "BedNorth"; case BedSlopeXEnum : return "BedSlopeX"; case BedSlopeYEnum : return "BedSlopeY"; @@ -712,6 +718,6 @@ case RheologyBbarAbsGradientEnum : return "RheologyBbarAbsGradient"; case SealevelEnum : return "Sealevel"; - case SealevelEustaticMaskEnum : return "SealevelEustaticMask"; - case SealevelEustaticOceanMaskEnum : return "SealevelEustaticOceanMask"; + case SealevelBarystaticMaskEnum : return "SealevelBarystaticMask"; + case SealevelBarystaticOceanMaskEnum : return "SealevelBarystaticOceanMask"; case SealevelNEsaEnum : return "SealevelNEsa"; case SealevelNEsaRateEnum : return "SealevelNEsaRate"; @@ -722,6 +728,5 @@ case BslrRateEnum : return "BslrRate"; case GmtslrEnum : return "Gmtslr"; - case SealevelRSLEustaticEnum : return "SealevelRSLEustatic"; - case SealevelRSLEustaticRateEnum : return "SealevelRSLEustaticRate"; + case SealevelRSLBarystaticEnum : return "SealevelRSLBarystatic"; case SealevelRSLRateEnum : return "SealevelRSLRate"; case SealevelUEastEsaEnum : return "SealevelUEastEsa"; @@ -729,14 +734,14 @@ case SealevelUEsaRateEnum : return "SealevelUEsaRate"; case SealevelUNorthEsaEnum : return "SealevelUNorthEsa"; - case SealevelriseCumDeltathicknessEnum : return "SealevelriseCumDeltathickness"; - case SealevelriseCumDeltathicknessOldEnum : return "SealevelriseCumDeltathicknessOld"; + case SealevelchangeCumDeltathicknessEnum : return "SealevelchangeCumDeltathickness"; + case SealevelchangeCumDeltathicknessOldEnum : return "SealevelchangeCumDeltathicknessOld"; case SurfaceloadOtherEnum : return "SurfaceloadOther"; case SurfaceloadIceThicknessChangeEnum : return "SurfaceloadIceThicknessChange"; case SurfaceloadWaterHeightChangeEnum : return "SurfaceloadWaterHeightChange"; - case SealevelriseIndicesEnum : return "SealevelriseIndices"; - case SealevelriseGEnum : return "SealevelriseG"; - case SealevelriseGUEnum : return "SealevelriseGU"; - case SealevelriseGEEnum : return "SealevelriseGE"; - case SealevelriseGNEnum : return "SealevelriseGN"; + case SealevelchangeIndicesEnum : return "SealevelchangeIndices"; + case SealevelchangeGEnum : return "SealevelchangeG"; + case SealevelchangeGUEnum : return "SealevelchangeGU"; + case SealevelchangeGEEnum : return "SealevelchangeGE"; + case SealevelchangeGNEnum : return "SealevelchangeGN"; case SedimentHeadEnum : return "SedimentHead"; case SedimentHeadOldEnum : return "SedimentHeadOld"; @@ -846,4 +851,9 @@ case SmbZTopEnum : return "SmbZTop"; case SmbZYEnum : return "SmbZY"; + case SolidearthExternalDisplacementEastRateEnum : return "SolidearthExternalDisplacementEastRate"; + case SolidearthExternalDisplacementNorthRateEnum : return "SolidearthExternalDisplacementNorthRate"; + case SolidearthExternalDisplacementUpRateEnum : return "SolidearthExternalDisplacementUpRate"; + case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate"; + case SolidearthExternalBarystaticSeaLevelRateEnum : return "SolidearthExternalBarystaticSeaLevelRate"; case StrainRateeffectiveEnum : return "StrainRateeffective"; case StrainRateparallelEnum : return "StrainRateparallel"; @@ -1305,6 +1315,6 @@ case SealevelNmotionEnum : return "SealevelNmotion"; case SealevelUmotionEnum : return "SealevelUmotion"; - case SealevelriseAnalysisEnum : return "SealevelriseAnalysis"; - case SealevelriseSolutionEnum : return "SealevelriseSolution"; + case SealevelchangeAnalysisEnum : return "SealevelchangeAnalysis"; + case SealevelchangeSolutionEnum : return "SealevelchangeSolution"; case SegEnum : return "Seg"; case SegInputEnum : return "SegInput"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25946) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25947) @@ -143,4 +143,8 @@ else if (strcmp(name,"DslModelid")==0) return DslModelidEnum; else if (strcmp(name,"DslNummodels")==0) return DslNummodelsEnum; + else if (strcmp(name,"SolidearthIsExternal")==0) return SolidearthIsExternalEnum; + else if (strcmp(name,"SolidearthExternalNature")==0) return SolidearthExternalNatureEnum; + else if (strcmp(name,"SolidearthExternalModelid")==0) return SolidearthExternalModelidEnum; + else if (strcmp(name,"SolidearthExternalNummodels")==0) return SolidearthExternalNummodelsEnum; else if (strcmp(name,"DslComputeFingerprints")==0) return DslComputeFingerprintsEnum; else if (strcmp(name,"EarthId")==0) return EarthIdEnum; @@ -256,12 +260,12 @@ else if (strcmp(name,"LoveFrequencies")==0) return LoveFrequenciesEnum; else if (strcmp(name,"LoveG0")==0) return LoveG0Enum; - else if (strcmp(name,"LoveKernels")==0) return LoveKernelsEnum; + else stage=3; + } + if(stage==3){ + if (strcmp(name,"LoveKernels")==0) return LoveKernelsEnum; else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum; else if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum; else if (strcmp(name,"LoveR0")==0) return LoveR0Enum; - else stage=3; - } - if(stage==3){ - if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum; + else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum; else if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum; else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum; @@ -346,5 +350,5 @@ else if (strcmp(name,"RotationalAngularVelocity")==0) return RotationalAngularVelocityEnum; else if (strcmp(name,"SolidearthSettingsElastic")==0) return SolidearthSettingsElasticEnum; - else if (strcmp(name,"SealevelriseGeometryDone")==0) return SealevelriseGeometryDoneEnum; + else if (strcmp(name,"SealevelchangeGeometryDone")==0) return SealevelchangeGeometryDoneEnum; else if (strcmp(name,"RotationalEquatorialMoi")==0) return RotationalEquatorialMoiEnum; else if (strcmp(name,"TidalLoveH")==0) return TidalLoveHEnum; @@ -355,11 +359,11 @@ else if (strcmp(name,"LoadLoveK")==0) return LoadLoveKEnum; else if (strcmp(name,"LoadLoveL")==0) return LoadLoveLEnum; - else if (strcmp(name,"SealevelriseGRigid")==0) return SealevelriseGRigidEnum; - else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum; + else if (strcmp(name,"SealevelchangeGRigid")==0) return SealevelchangeGRigidEnum; + else if (strcmp(name,"SealevelchangeGElastic")==0) return SealevelchangeGElasticEnum; else if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum; else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum; else if (strcmp(name,"SolidearthSettingsGlfraction")==0) return SolidearthSettingsGlfractionEnum; else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum; - else if (strcmp(name,"SealevelriseHElastic")==0) return SealevelriseHElasticEnum; + else if (strcmp(name,"SealevelchangeHElastic")==0) return SealevelchangeHElasticEnum; else if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum; else if (strcmp(name,"SolidearthSettingsMaxiter")==0) return SolidearthSettingsMaxiterEnum; @@ -367,10 +371,10 @@ else if (strcmp(name,"RotationalPolarMoi")==0) return RotationalPolarMoiEnum; else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum; - else if (strcmp(name,"SealevelriseRequestedOutputs")==0) return SealevelriseRequestedOutputsEnum; + else if (strcmp(name,"SealevelchangeRequestedOutputs")==0) return SealevelchangeRequestedOutputsEnum; else if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum; else if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum; - else if (strcmp(name,"SealevelriseRunCount")==0) return SealevelriseRunCountEnum; - else if (strcmp(name,"SealevelriseTransitions")==0) return SealevelriseTransitionsEnum; - else if (strcmp(name,"SealevelriseUElastic")==0) return SealevelriseUElasticEnum; + else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum; + else if (strcmp(name,"SealevelchangeTransitions")==0) return SealevelchangeTransitionsEnum; + else if (strcmp(name,"SealevelchangeUElastic")==0) return SealevelchangeUElasticEnum; else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum; else if (strcmp(name,"SettingsNumResultsOnNodes")==0) return SettingsNumResultsOnNodesEnum; @@ -379,12 +383,12 @@ else if (strcmp(name,"SettingsResultsOnNodes")==0) return SettingsResultsOnNodesEnum; else if (strcmp(name,"SettingsSbCouplingFrequency")==0) return SettingsSbCouplingFrequencyEnum; - else if (strcmp(name,"SettingsSolverResidueThreshold")==0) return SettingsSolverResidueThresholdEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"SettingsSolverResidueThreshold")==0) return SettingsSolverResidueThresholdEnum; else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum; else if (strcmp(name,"SmbAIce")==0) return SmbAIceEnum; else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum; + else if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum; else if (strcmp(name,"SmbAccualti")==0) return SmbAccualtiEnum; else if (strcmp(name,"SmbAccugrad")==0) return SmbAccugradEnum; @@ -502,12 +506,12 @@ else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum; else if (strcmp(name,"Velocity")==0) return VelocityEnum; - else if (strcmp(name,"WorldComm")==0) return WorldCommEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"WorldComm")==0) return WorldCommEnum; else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum; else if (strcmp(name,"InputsSTART")==0) return InputsSTARTEnum; else if (strcmp(name,"Adjoint")==0) return AdjointEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"Adjointp")==0) return AdjointpEnum; + else if (strcmp(name,"Adjointp")==0) return AdjointpEnum; else if (strcmp(name,"Adjointx")==0) return AdjointxEnum; else if (strcmp(name,"Adjointy")==0) return AdjointyEnum; @@ -546,4 +550,6 @@ else if (strcmp(name,"BaseSlopeY")==0) return BaseSlopeYEnum; else if (strcmp(name,"Bed")==0) return BedEnum; + else if (strcmp(name,"BedEast")==0) return BedEastEnum; + else if (strcmp(name,"BedNorth")==0) return BedNorthEnum; else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum; else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum; @@ -623,5 +629,8 @@ else if (strcmp(name,"FrictionCoefficientcoulomb")==0) return FrictionCoefficientcoulombEnum; else if (strcmp(name,"FrictionEffectivePressure")==0) return FrictionEffectivePressureEnum; - else if (strcmp(name,"FrictionM")==0) return FrictionMEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"FrictionM")==0) return FrictionMEnum; else if (strcmp(name,"FrictionP")==0) return FrictionPEnum; else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum; @@ -629,8 +638,5 @@ else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum; else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; + else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; else if (strcmp(name,"Frictionf")==0) return FrictionfEnum; else if (strcmp(name,"FrontalForcingsBasinId")==0) return FrontalForcingsBasinIdEnum; @@ -727,6 +733,6 @@ else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; else if (strcmp(name,"Sealevel")==0) return SealevelEnum; - else if (strcmp(name,"SealevelEustaticMask")==0) return SealevelEustaticMaskEnum; - else if (strcmp(name,"SealevelEustaticOceanMask")==0) return SealevelEustaticOceanMaskEnum; + else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum; + else if (strcmp(name,"SealevelBarystaticOceanMask")==0) return SealevelBarystaticOceanMaskEnum; else if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum; else if (strcmp(name,"SealevelNEsaRate")==0) return SealevelNEsaRateEnum; @@ -737,6 +743,5 @@ else if (strcmp(name,"BslrRate")==0) return BslrRateEnum; else if (strcmp(name,"Gmtslr")==0) return GmtslrEnum; - else if (strcmp(name,"SealevelRSLEustatic")==0) return SealevelRSLEustaticEnum; - else if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum; + else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum; else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum; else if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum; @@ -744,17 +749,17 @@ else if (strcmp(name,"SealevelUEsaRate")==0) return SealevelUEsaRateEnum; else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum; - else if (strcmp(name,"SealevelriseCumDeltathickness")==0) return SealevelriseCumDeltathicknessEnum; - else if (strcmp(name,"SealevelriseCumDeltathicknessOld")==0) return SealevelriseCumDeltathicknessOldEnum; + else if (strcmp(name,"SealevelchangeCumDeltathickness")==0) return SealevelchangeCumDeltathicknessEnum; + else if (strcmp(name,"SealevelchangeCumDeltathicknessOld")==0) return SealevelchangeCumDeltathicknessOldEnum; else if (strcmp(name,"SurfaceloadOther")==0) return SurfaceloadOtherEnum; - else if (strcmp(name,"SurfaceloadIceThicknessChange")==0) return SurfaceloadIceThicknessChangeEnum; - else if (strcmp(name,"SurfaceloadWaterHeightChange")==0) return SurfaceloadWaterHeightChangeEnum; - else if (strcmp(name,"SealevelriseIndices")==0) return SealevelriseIndicesEnum; - else if (strcmp(name,"SealevelriseG")==0) return SealevelriseGEnum; - else if (strcmp(name,"SealevelriseGU")==0) return SealevelriseGUEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SealevelriseGE")==0) return SealevelriseGEEnum; - else if (strcmp(name,"SealevelriseGN")==0) return SealevelriseGNEnum; + if (strcmp(name,"SurfaceloadIceThicknessChange")==0) return SurfaceloadIceThicknessChangeEnum; + else if (strcmp(name,"SurfaceloadWaterHeightChange")==0) return SurfaceloadWaterHeightChangeEnum; + else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum; + else if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum; + else if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum; + else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum; + else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum; else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; @@ -864,6 +869,14 @@ else if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum; else if (strcmp(name,"SmbZY")==0) return SmbZYEnum; + else if (strcmp(name,"SolidearthExternalDisplacementEastRate")==0) return SolidearthExternalDisplacementEastRateEnum; + else if (strcmp(name,"SolidearthExternalDisplacementNorthRate")==0) return SolidearthExternalDisplacementNorthRateEnum; + else if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum; + else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum; + else if (strcmp(name,"SolidearthExternalBarystaticSeaLevelRate")==0) return SolidearthExternalBarystaticSeaLevelRateEnum; else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum; - else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum; else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; @@ -875,8 +888,5 @@ else if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum; else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; + else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; @@ -988,5 +998,8 @@ else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum; else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum; - else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum; else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum; else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum; @@ -998,8 +1011,5 @@ else if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum; else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum; + else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum; else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum; else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum; @@ -1111,5 +1121,8 @@ else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum; else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum; - else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; else if (strcmp(name,"Element")==0) return ElementEnum; @@ -1121,8 +1134,5 @@ else if (strcmp(name,"EsaTransitions")==0) return EsaTransitionsEnum; else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum; + else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum; else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum; @@ -1234,5 +1244,8 @@ else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum; else if (strcmp(name,"Materials")==0) return MaterialsEnum; - else if (strcmp(name,"Matestar")==0) return MatestarEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"Matestar")==0) return MatestarEnum; else if (strcmp(name,"Matice")==0) return MaticeEnum; else if (strcmp(name,"Matlitho")==0) return MatlithoEnum; @@ -1244,8 +1257,5 @@ else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum; else if (strcmp(name,"MaxVel")==0) return MaxVelEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"MaxVx")==0) return MaxVxEnum; + else if (strcmp(name,"MaxVx")==0) return MaxVxEnum; else if (strcmp(name,"MaxVy")==0) return MaxVyEnum; else if (strcmp(name,"MaxVz")==0) return MaxVzEnum; @@ -1335,6 +1345,6 @@ else if (strcmp(name,"SealevelNmotion")==0) return SealevelNmotionEnum; else if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum; - else if (strcmp(name,"SealevelriseAnalysis")==0) return SealevelriseAnalysisEnum; - else if (strcmp(name,"SealevelriseSolution")==0) return SealevelriseSolutionEnum; + else if (strcmp(name,"SealevelchangeAnalysis")==0) return SealevelchangeAnalysisEnum; + else if (strcmp(name,"SealevelchangeSolution")==0) return SealevelchangeSolutionEnum; else if (strcmp(name,"Seg")==0) return SegEnum; else if (strcmp(name,"SegInput")==0) return SegInputEnum; @@ -1357,5 +1367,8 @@ else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum; else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum; - else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum; else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum; else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum; @@ -1367,8 +1380,5 @@ else if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum; else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; + else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; Index: /issm/trunk-jpl/src/m/classes/additionalsolidearthsolution.m =================================================================== --- /issm/trunk-jpl/src/m/classes/additionalsolidearthsolution.m (revision 25947) +++ /issm/trunk-jpl/src/m/classes/additionalsolidearthsolution.m (revision 25947) @@ -0,0 +1,43 @@ +%ADDITIONALSOLIDEARTHSOLUTION class definition +% +% Usage: +% addsol=additionalsolidearthsolution(); + +classdef additionalsolidearthsolution < solidearthsolution + properties (SetAccess=public) + end + methods + function self = additionalsolidearthsolution(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + otherwise + error('constructor not supported'); + end + end % }}} + function self = setdefaultparameters(self) % {{{ + self.setdefaultparameters@solidearthsolution(); + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + if ~ismember('SealevelriseAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.solidearth.settings.isgrd==0), + error('additionalsolidearthsolution checkconsistency error message: need to run GRD solution if you are supplying a GRD additional pattern solution'); + end + self.checkconsistency@solidearthsolution(md,solution,analyses); + + end % }}} + function disp(self) % {{{ + disp(sprintf(' external: additionalsolidearth solution:')); + self.disp@solidearthsolution(); + end % }}} + function marshall(self,prefix,md,fid) % {{{ + WriteData(fid,prefix,'object',self,'data',1,'name','md.solidearth.external.nature','format','Integer'); %code 1 for additionalsolidearthsolution class + self.marshall@solidearthsolution(prefix,md,fid); + end % }}} + function savemodeljs(self,fid,modelname) % {{{ + self.savemodeljs@solidearthsolution(fid,modelname); + end % }}} + function self = extrude(self,md) % {{{ + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/mmeadditionalsolidearthsolution.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mmeadditionalsolidearthsolution.m (revision 25947) +++ /issm/trunk-jpl/src/m/classes/mmeadditionalsolidearthsolution.m (revision 25947) @@ -0,0 +1,86 @@ +%MMEADDITIONALSOLIDEARTHSOLUTION class definition +% +% Usage: +% addsol=mmeadditionalsolidearthsolution(); where the additional solid earth solutions % are based on a multi-model ensemble (ex: Caron et al 2017 statistics) + +classdef mmeadditionalsolidearthsolution < additionalsolidearthsolution + properties (SetAccess=public) + modelid; %index into the multi-model ensemble, each ensemble variable being defined + %in the father class. + end + methods + function self = mmeadditionalsolidearthsolution(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + otherwise + error('constructor not supported'); + end + end % }}} + function self = setdefaultparameters(self) % {{{ + self.setdefaultparameters@additionalsolidearthsolution(); + self.modelid=0; + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + if ~ismember('SealevelriseAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.solidearth.settings.isgrd==0), + error('mmeadditionalsolidearthsolution checkconsistency error message: need to run GRD solution if you are supplying a GRD additional pattern solution'); + end + + seast=length(self.displacementeast); + snorth=length(self.displacementnorth); + sup=length(self.displacementup); + sgeoid=length(self.geoid); + sbarystaticsealevel=length(self.barystaticsealevel); + + if (seast-snorth)~=0, + error('mmeadditionalsolidearthsolution checkconsistency error message: displacementeast and displacementnorth should be the same size'); + end + + if (seast-sup)~=0, + error('mmeadditionalsolidearthsolution checkconsistency error message: displacementeast and displacementup should be the same size'); + end + + if (seast-sgeoid)~=0, + error('mmeadditionalsolidearthsolution checkconsistency error message: displacementeast and geoid should be the same size'); + end + + if (seast-sbarystaticsealevel)~=0, + error('mmeadditionalsolidearthsolution checkconsistency error message: displacementeast and barystaticsealevel should be the same size'); + end + + md = checkfield(md,'field',self.modelid,'NaN',1,'Inf',1,'>=',1,'<=',length(self.displacementeast)); + + for i=1:seast, + md = checkfield(md,'field',self.displacementeast{i},'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'field',self.displacementup{i},'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'field',self.displacementnorth{i},'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'field',self.geoid{i},'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'field',self.barystaticsealevel{i},'NaN',1,'Inf',1,'timeseries',1); + end + + + end % }}} + function disp(self) % {{{ + disp(sprintf(' external: mmeadditionalsolidearth solution:')); + self.disp@solidearthsolution(); + fielddisplay(self,'modelid','index into the multi-model ensemble, determines which field will be used.'); + end % }}} + function marshall(self,prefix,md,fid) % {{{ + WriteData(fid,prefix,'object',self,'data',3,'name','md.solidearth.external.nature','format','Integer'); %code 3 for mmeadditionalsolidearthsolution class + WriteData(fid,prefix,'name','md.solidearth.external.nummodels','data',length(self.displacementeast),'format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','modelid','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','displacementeast','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); + WriteData(fid,prefix,'object',self,'fieldname','displacementup','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); + WriteData(fid,prefix,'object',self,'fieldname','displacementnorth','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); + WriteData(fid,prefix,'object',self,'fieldname','geoid','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); + WriteData(fid,prefix,'object',self,'fieldname','barystaticsealevel','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); + + end % }}} + function savemodeljs(self,fid,modelname) % {{{ + error('mmeadditionalsolidearthsolution error message: not implemented yet'); + end % }}} + function self = extrude(self,md) % {{{ + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/mmeofflinesolidearthsolution.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mmeofflinesolidearthsolution.m (revision 25947) +++ /issm/trunk-jpl/src/m/classes/mmeofflinesolidearthsolution.m (revision 25947) @@ -0,0 +1,86 @@ +%MMEOFFLINESOLIDEARTHSOLUTION class definition +% +% Usage: +% addsol=mmeofflinesolidearthsolution(); where the offline solid earth solutions % are based on a multi-model ensemble (ex: Caron et al 2017 statistics) + +classdef mmeofflinesolidearthsolution < offlinesolidearthsolution + properties (SetAccess=public) + modelid; %index into the multi-model ensemble, each ensemble variable being defined + %in the father class. + end + methods + function self = mmeofflinesolidearthsolution(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + otherwise + error('constructor not supported'); + end + end % }}} + function self = setdefaultparameters(self) % {{{ + self.setdefaultparameters@offlinesolidearthsolution(); + self.modelid=0; + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + if ~ismember('SealevelriseAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.solidearth.settings.isgrd==1), + error('mmeofflinesolidearthsolution checkconsistency error message: trying to run GRD patterns while supplying an offline solution for those patterns!'); + end + + seast=length(self.displacementeast); + snorth=length(self.displacementnorth); + sup=length(self.displacementup); + sgeoid=length(self.geoid); + sbarystaticsealevel=length(self.barystaticsealevel); + + if (seast-snorth)~=0, + error('mmeofflinesolidearthsolution checkconsistency error message: displacementeast and displacementnorth should be the same size'); + end + + if (seast-sup)~=0, + error('mmeofflinesolidearthsolution checkconsistency error message: displacementeast and displacementup should be the same size'); + end + + if (seast-sgeoid)~=0, + error('mmeofflinesolidearthsolution checkconsistency error message: displacementeast and geoid should be the same size'); + end + + if (seast-sbarystaticsealevel)~=0, + error('mmeofflinesolidearthsolution checkconsistency error message: displacementeast and barystaticsealevel should be the same size'); + end + + md = checkfield(md,'field',self.modelid,'NaN',1,'Inf',1,'>=',1,'<=',length(self.displacementeast)); + + for i=1:seast, + md = checkfield(md,'field',self.displacementeast{i},'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'field',self.displacementup{i},'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'field',self.displacementnorth{i},'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'field',self.geoid{i},'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'field',self.barystaticsealevel{i},'NaN',1,'Inf',1,'timeseries',1); + end + + + end % }}} + function disp(self) % {{{ + disp(sprintf(' external: mmeofflinesolidearth solution:')); + self.disp@solidearthsolution(); + fielddisplay(self,'modelid','index into the multi-model ensemble, determines which field will be used.'); + end % }}} + function marshall(self,prefix,md,fid) % {{{ + WriteData(fid,prefix,'object',self,'data',4,'name','md.solidearth.external.nature','format','Integer'); %code 4 for mmeofflinesolidearthsolution class + WriteData(fid,prefix,'object',self,'fieldname','modelid','format','Double'); + WriteData(fid,prefix,'name','md.solidearth.external.nummodels','data',length(self.displacementeast),'format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','displacementeast','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); + WriteData(fid,prefix,'object',self,'fieldname','displacementup','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); + WriteData(fid,prefix,'object',self,'fieldname','displacementnorth','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); + WriteData(fid,prefix,'object',self,'fieldname','geoid','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); + WriteData(fid,prefix,'object',self,'fieldname','barystaticsealevel','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); + + end % }}} + function savemodeljs(self,fid,modelname) % {{{ + error('mmeofflinesolidearthsolution error message: not implemented yet'); + end % }}} + function self = extrude(self,md) % {{{ + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/offlinesolidearthsolution.m =================================================================== --- /issm/trunk-jpl/src/m/classes/offlinesolidearthsolution.m (revision 25947) +++ /issm/trunk-jpl/src/m/classes/offlinesolidearthsolution.m (revision 25947) @@ -0,0 +1,44 @@ +%OFFLINESOLIDEARTHSOLUTION class definition +% +% Usage: +% addsol=offlinesolidearthsolution(); + +classdef offlinesolidearthsolution < solidearthsolution + properties (SetAccess=public) + end + methods + function self = offlinesolidearthsolution(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + otherwise + error('constructor not supported'); + end + end % }}} + function self = setdefaultparameters(self) % {{{ + self.setdefaultparameters@solidearthsolution(); + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + if ~ismember('SealevelriseAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.solidearth.settings.isgrd==1), + error('offlinesolidearthsolution checkconsistency error message: trying to run GRD patterns while supplying an offline solution for those patterns!'); + end + self.checkconsistency@solidearthsolution(md,solution,analyses); + + end % }}} + function disp(self) % {{{ + disp(sprintf(' external: offlinesolidearth solution:')); + self.disp@solidearthsolution(); + end % }}} + function marshall(self,prefix,md,fid) % {{{ + disp(sprintf(' external (offlinesolidearthsolution) solution:')); + WriteData(fid,prefix,'object',self,'data',2,'name','md.solidearth.external.nature','format','Integer'); %code 1 for offlinesolidearthsolution class + self.marshall@solidearthsolution(prefix,md,fid); + end % }}} + function savemodeljs(self,fid,modelname) % {{{ + self.savemodeljs@solidearthsolution(fid,modelname); + end % }}} + function self = extrude(self,md) % {{{ + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/solidearth.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearth.m (revision 25946) +++ /issm/trunk-jpl/src/m/classes/solidearth.m (revision 25947) @@ -41,5 +41,5 @@ self.partitionhydro=[]; - %no external solutions by defalt: + %no external solutions by default: self.external=[]; @@ -120,5 +120,8 @@ self.rotational.marshall(prefix,md,fid); if ~isempty(self.external), + WriteData(fid,prefix,'data',1,'format','Integer','name','md.solidearth.isexternal'); self.external.marshall(prefix,md,fid); + else + WriteData(fid,prefix,'data',0,'format','Integer','name','md.solidearth.isexternal'); end Index: /issm/trunk-jpl/src/m/classes/solidearthsolution.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsolution.m (revision 25946) +++ /issm/trunk-jpl/src/m/classes/solidearthsolution.m (revision 25947) @@ -40,18 +40,30 @@ end % }}} function disp(self) % {{{ - fielddisplay(self,'displacementeast','solid-Earth Eastwards bedrock motion (m/yr)'); - fielddisplay(self,'displacementnorth','solid-Earth Northwards bedrock motion (m/yr)'); - fielddisplay(self,'displacementup','solid-Earth bedrock uplift rate (m/yr)'); - fielddisplay(self,'geoid','solid-Earth geoid motion (m/yr)'); - fielddisplay(self,'barystaticsealevel','barystatic absolute sea-level rate (m/yr)'); + disp(sprintf('%s',' units for time series is (yr)')); + fielddisplay(self,'displacementeast','solid-Earth Eastwards bedrock displacement series (m)'); + fielddisplay(self,'displacementnorth','solid-Earth Northwards bedrock displacement time series (m)'); + fielddisplay(self,'displacementup','solid-Earth bedrock uplift time series (m)'); + fielddisplay(self,'geoid','solid-Earth geoid time series (m)'); + fielddisplay(self,'barystaticsealevel','barystatic absolute sea-level time series (m)'); end % }}} function marshall(self,prefix,md,fid) % {{{ yts=md.constants.yts; - WriteData(fid,prefix,'object',self,'class','solidearthsolution','fieldname','displacementeast','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'class','solidearthsolution','fieldname','displacementnorth','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'class','solidearthsolution','fieldname','displacementup','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'class','solidearthsolution','fieldname','geoid','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'class','solidearthsolution','fieldname','barystaticsealevel','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + + %transform our time series into time series rates + time=self.displacementeast(end,:); + dt=diff(time,1,2); + + displacementeast_rate=diff(displacementeast(1:end-1,:),1,2)./dt; + displacementnorth_rate=diff(displacementnorth(1:end-1,:),1,2)./dt; + displacementup_rate=diff(displacementup(1:end-1,:),1,2)./dt; + geoid_rate=diff(geoid(1:end-1,:),1,2)./dt; + barystaticsealevel_rate=diff(barystaticsealevel(1:end-1,:),1,2)./dt; + + WriteData(fid,prefix,'name','md.solidearth.external.displacementeast','data',displacementeast_rate,'format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'name','md.solidearth.external.displacementup','data',displacementup_rate,'format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'name','md.solidearth.external.displacementnorth','data',displacementnorth_rate,'format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'name','md.solidearth.external.geoid','data',geoid_rate,'format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'name','md.solidearth.external.barystaticsealevel','data',barystaticsealevel_rate,'format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); end % }}}