Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 25763) @@ -32,5 +32,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum); Index: /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 25763) @@ -99,5 +99,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum); Index: /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp (revision 25763) @@ -41,5 +41,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum); Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 25763) @@ -136,5 +136,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.thickness",ThicknessEnum); iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 25763) @@ -74,5 +74,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum,0.); Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 25763) @@ -77,5 +77,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum); Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 25763) @@ -129,5 +129,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); /*Needed for friction*/ + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); /*Needed for friction*/ iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum); iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum); Index: /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 25763) @@ -58,5 +58,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); if(iomodel->domaintype!=Domain2DhorizontalEnum){ iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); Index: /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 25763) @@ -42,5 +42,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); if(iomodel->domaintype!=Domain2DhorizontalEnum & iomodel->domaintype!=Domain3DsurfaceEnum){ Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 25763) @@ -148,5 +148,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum); Index: /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp (revision 25763) @@ -56,5 +56,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); if(iomodel->domaintype!=Domain2DhorizontalEnum){ Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 25763) @@ -40,5 +40,5 @@ iomodel->FetchData(&geodetic,"md.solidearth.settings.computesealevelchange"); iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.icethicknesschange",SurfaceloadIceThicknessChangeEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + 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); @@ -175,5 +175,7 @@ IssmDouble planetradius=0; IssmDouble planetarea=0; + bool rigid=false; bool elastic=false; + bool rotation=false; int numoutputs; @@ -207,4 +209,5 @@ 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)); @@ -260,7 +263,16 @@ } /*}}}*/ /*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){ - /*love numbers: */ iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h"); iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k"); @@ -269,5 +281,4 @@ iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk"); iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl"); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1)); @@ -278,44 +289,63 @@ parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1)); - /*compute elastic green function for a range of angles*/ - iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); - M=reCast(180./degacc+1.); - // AD performance is sensitive to calls to ensurecontiguous. // // Providing "t" will cause ensurecontiguous to be called. #ifdef _HAVE_AD_ - G_rigid=xNew(M,"t"); G_elastic=xNew(M,"t"); U_elastic=xNew(M,"t"); H_elastic=xNew(M,"t"); #else - G_rigid=xNew(M); 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()); - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. + } + if(elastic){ #ifdef _HAVE_AD_ G_elastic_local=xNew(m,"t"); - G_rigid_local=xNew(m,"t"); U_elastic_local=xNew(m,"t"); H_elastic_local=xNew(m,"t"); #else G_elastic_local=xNew(m); - G_rigid_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]; @@ -356,4 +386,6 @@ } } + } + if(rigid){ /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/ @@ -369,22 +401,22 @@ /*All gather:*/ ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - 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()); - + 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]; - parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); - G_elastic[0]=G_elastic[1]; - parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M)); - U_elastic[0]=U_elastic[1]; - parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M)); - H_elastic[0]=H_elastic[1]; - parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M)); + 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*/ @@ -397,19 +429,21 @@ /*free resources: */ - xDelete(love_h); - xDelete(love_k); - xDelete(love_l); - xDelete(love_th); - xDelete(love_tk); - xDelete(love_tl); xDelete(G_rigid); xDelete(G_rigid_local); - xDelete(G_elastic); - xDelete(G_elastic_local); - xDelete(U_elastic); - xDelete(U_elastic_local); - xDelete(H_elastic); - xDelete(H_elastic_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: */ @@ -436,5 +470,4 @@ /*}}}*/ - /*}}}*/ }/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 25763) @@ -768,5 +768,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum); Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 25763) @@ -105,5 +105,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); if(iomodel->domaintype!=Domain2DhorizontalEnum){ Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 25763) @@ -133,5 +133,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); - iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25763) @@ -5463,5 +5463,5 @@ this->inputs->SetDoubleInput(AreaEnum,this->lid,area); this->AddInput(SealevelAreaEnum,&area,P0Enum); - if(!computerigid & !computeelastic)return; + if(!computerigid)return; /*recover precomputed green function kernels:*/ @@ -5469,23 +5469,27 @@ parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M); - parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M); - - parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M); - - parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); + if(computeelastic){ + parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M); + + parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M); + + parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); + } /*allocate indices:*/ indices=xNew(gsize); G=xNewZeroInit(gsize); - GU=xNewZeroInit(gsize); - if(horiz){ - GN=xNewZeroInit(gsize); - GE=xNewZeroInit(gsize); - } - - /* Where is the centroid of this element:{{{*/ + if(computeelastic){ + GU=xNewZeroInit(gsize); + if(horiz){ + GN=xNewZeroInit(gsize); + GE=xNewZeroInit(gsize); + } + } + + /* Where is the centroid of this element: /*retrieve coordinates: */ @@ -5515,7 +5519,6 @@ /*Rigid earth gravitational perturbation: */ - if(computerigid){ - G[i]+=G_rigid_precomputed[index]; - } + G[i]+=G_rigid_precomputed[index]; + if(computeelastic){ G[i]+=G_elastic_precomputed[index]; @@ -5524,20 +5527,22 @@ /*Elastic components:*/ - GU[i] = constant * U_elastic_precomputed[index]; - if(horiz){ - /*Compute azimuths, both north and east components: */ - x = xx[i]; y = yy[i]; z = zz[i]; - if(latitude[i]==90){ - x=1e-12; y=1e-12; - } - if(latitude[i]==-90){ - x=1e-12; y=1e-12; - } - dx = x_element-x; dy = y_element-y; dz = z_element-z; - N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); - E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); - - GN[i] = constant*H_elastic_precomputed[index]*N_azim; - GE[i] = constant*H_elastic_precomputed[index]*E_azim; + if(computeelastic){ + GU[i] = constant * U_elastic_precomputed[index]; + if(horiz){ + /*Compute azimuths, both north and east components: */ + x = xx[i]; y = yy[i]; z = zz[i]; + if(latitude[i]==90){ + x=1e-12; y=1e-12; + } + if(latitude[i]==-90){ + x=1e-12; y=1e-12; + } + dx = x_element-x; dy = y_element-y; dz = z_element-z; + N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + + GN[i] = constant*H_elastic_precomputed[index]*N_azim; + GE[i] = constant*H_elastic_precomputed[index]*E_azim; + } } } @@ -5546,8 +5551,10 @@ this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); - this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize); - if(horiz){ - this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize); - this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize); + if(computeelastic){ + this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize); + if(horiz){ + this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize); + this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize); + } } @@ -5556,16 +5563,20 @@ delete indices; delete G; - delete GU; - if(horiz){ - delete GN; - delete GE; + if(computeelastic){ + delete GU; + if(horiz){ + delete GN; + delete GE; + } } #else xDelete(indices); xDelete(G); - xDelete(GU); - if(horiz){ - xDelete(GN); - xDelete(GE); + if(computeelastic){ + xDelete(GU); + if(horiz){ + xDelete(GN); + xDelete(GE); + } } #endif Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 25763) @@ -24,4 +24,5 @@ bool isslr=0; bool isgia=0; + int grd=0; int solution_type; @@ -30,4 +31,5 @@ femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/ @@ -59,5 +61,5 @@ /*Run geodetic:*/ - grd_core(femmodel); + if(grd)grd_core(femmodel); /*Run steric core for sure:*/ @@ -109,14 +111,11 @@ int frequency,count; int horiz; - int geodetic=0; IssmDouble dt; IssmDouble oceanarea; int bp_compute_fingerprints=0; - /*Should we even be here?:*/ - femmodel->parameters->FindParam(&geodetic,SolidearthSettingsComputesealevelchangeEnum); if(!geodetic)return; - + /*Verbose: */ - if(VerboseSolution()) _printf0_(" computing geodetic sea level rise\n"); + if(VerboseSolution()) _printf0_(" computing GRD sea level patterns\n"); /*retrieve more parameters:*/ @@ -178,8 +177,8 @@ RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,oceanarea); - /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */ + /*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) geodetic signatures:*/ + /*compute viscosus (GIA) signatures:*/ sealevelrise_core_viscous(&U_gia,&N_gia,femmodel,RSLg); @@ -269,5 +268,5 @@ int solution_type; IssmDouble dt; - int geodetic=0; + int grd=0; int step; IssmDouble time; @@ -281,5 +280,5 @@ femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); - femmodel->parameters->FindParam(&geodetic,SolidearthSettingsComputesealevelchangeEnum); + femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); femmodel->parameters->FindParam(&cumgmtslr,CumGmtslrEnum); femmodel->parameters->FindParam(&cumbslr,CumBslrEnum); @@ -301,5 +300,5 @@ GetStericRate(&steric_rate_g,femmodel); GetDynamicRate(&dynamic_rate_g,femmodel); - if(geodetic){ + if(grd){ GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum); GetVectorFromInputsx(&U_gia_rate,femmodel,UGiaRateEnum,VertexSIdEnum); @@ -318,5 +317,5 @@ /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate) * dt + steric_rate + dynamic_rate dt*/ - if(geodetic){ + if(grd){ SL->AXPY(N_gia_rate,dt); SL->AXPY(N_esa_rate,dt); @@ -326,5 +325,5 @@ /*compute new bedrock position: */ - if(geodetic){ + if(grd){ bedrock->AXPY(U_esa_rate,dt); bedrock->AXPY(U_gia_rate,dt); @@ -340,5 +339,5 @@ delete steric_rate_g; delete dynamic_rate_g; - if(geodetic){ + if(grd){ delete U_esa_rate; delete U_gia_rate; Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25762) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25763) @@ -351,4 +351,5 @@ syn keyword cConstant SealevelriseGElasticEnum syn keyword cConstant SolidearthSettingsComputesealevelchangeEnum +syn keyword cConstant SolidearthSettingsGRDEnum syn keyword cConstant SolidearthSettingsGlfractionEnum syn keyword cConstant SolidearthSettingsRunFrequencyEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25762) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25763) @@ -345,4 +345,5 @@ SealevelriseGElasticEnum, SolidearthSettingsComputesealevelchangeEnum, + SolidearthSettingsGRDEnum, SolidearthSettingsGlfractionEnum, SolidearthSettingsRunFrequencyEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25763) @@ -353,4 +353,5 @@ case SealevelriseGElasticEnum : return "SealevelriseGElastic"; case SolidearthSettingsComputesealevelchangeEnum : return "SolidearthSettingsComputesealevelchange"; + case SolidearthSettingsGRDEnum : return "SolidearthSettingsGRD"; case SolidearthSettingsGlfractionEnum : return "SolidearthSettingsGlfraction"; case SolidearthSettingsRunFrequencyEnum : return "SolidearthSettingsRunFrequency"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25762) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25763) @@ -359,4 +359,5 @@ else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum; 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; @@ -382,9 +383,9 @@ 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; + if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum; + else if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum; else if (strcmp(name,"SmbAccualti")==0) return SmbAccualtiEnum; else if (strcmp(name,"SmbAccugrad")==0) return SmbAccugradEnum; @@ -505,9 +506,9 @@ 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; + if (strcmp(name,"Adjoint")==0) return AdjointEnum; + else if (strcmp(name,"Adjointp")==0) return AdjointpEnum; else if (strcmp(name,"Adjointx")==0) return AdjointxEnum; else if (strcmp(name,"Adjointy")==0) return AdjointyEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum; else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum; - else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"Frictionf")==0) return FrictionfEnum; + if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; + else if (strcmp(name,"Frictionf")==0) return FrictionfEnum; else if (strcmp(name,"FrontalForcingsBasinId")==0) return FrontalForcingsBasinIdEnum; else if (strcmp(name,"FrontalForcingsSubglacialDischarge")==0) return FrontalForcingsSubglacialDischargeEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"SealevelriseGE")==0) return SealevelriseGEEnum; else if (strcmp(name,"SealevelriseGN")==0) return SealevelriseGNEnum; - else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; + if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; + else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; else if (strcmp(name,"SedimentHeadSubstep")==0) return SedimentHeadSubstepEnum; else if (strcmp(name,"SedimentHeadTransient")==0) return SedimentHeadTransientEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum; else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum; - else if (strcmp(name,"Area")==0) return AreaEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"SealevelArea")==0) return SealevelAreaEnum; + if (strcmp(name,"Area")==0) return AreaEnum; + else if (strcmp(name,"SealevelArea")==0) return SealevelAreaEnum; else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum; else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; @@ -997,9 +998,9 @@ else if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum; else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum; - else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum; else stage=9; } if(stage==9){ - if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum; + if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum; + else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum; else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum; else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum; else if (strcmp(name,"FemModel")==0) return FemModelEnum; - else if (strcmp(name,"FileParam")==0) return FileParamEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum; + if (strcmp(name,"FileParam")==0) return FileParamEnum; + else if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum; else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum; else if (strcmp(name,"FloatingAreaScaled")==0) return FloatingAreaScaledEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"MeshX")==0) return MeshXEnum; else if (strcmp(name,"MeshY")==0) return MeshYEnum; - else if (strcmp(name,"MinVel")==0) return MinVelEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"MinVx")==0) return MinVxEnum; + if (strcmp(name,"MinVel")==0) return MinVelEnum; + else if (strcmp(name,"MinVx")==0) return MinVxEnum; else if (strcmp(name,"MinVy")==0) return MinVyEnum; else if (strcmp(name,"MinVz")==0) return MinVzEnum; @@ -1366,9 +1367,9 @@ else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum; else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum; - else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum; + if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum; + else if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum; else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum; else if (strcmp(name,"TotalFloatingBmbScaled")==0) return TotalFloatingBmbScaledEnum; Index: /issm/trunk-jpl/src/m/classes/solidearth.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearth.m (revision 25762) +++ /issm/trunk-jpl/src/m/classes/solidearth.m (revision 25763) @@ -6,6 +6,7 @@ classdef solidearth properties (SetAccess=public) - sealevel = NaN; + initialsealevel = NaN; settings = solidearthsettings(); + external = []; surfaceload = surfaceload(); lovenumbers = lovenumbers(); @@ -38,4 +39,7 @@ self.partitionhydro=[]; + %no external solutions by defalt: + self.external=[]; + %earth radius self.planetradius= planetradius('earth'); @@ -48,5 +52,5 @@ end - md = checkfield(md,'fieldname','solidearth.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','solidearth.initialsealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','solidearth.requested_outputs','stringrow',1); @@ -55,4 +59,10 @@ self.lovenumbers.checkconsistency(md,solution,analyses); self.rotational.checkconsistency(md,solution,analyses); + if ~isempty(self.external), + if ~isa(self.external,'solidearthsolution'), + error('solidearth consistency check: external field should be a solidearthsolution'); + end + self.external.checkconsistency(md,solution,analyses); + end @@ -64,5 +74,5 @@ disp(sprintf(' solidearth inputs, forcings and settings:')); - fielddisplay(self,'sealevel','current sea level (prior to computation) [m]'); + fielddisplay(self,'initialsealevel','sea level at the start of computation) [m]'); fielddisplay(self,'planetradius','planet radius [m]'); fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'); @@ -70,13 +80,17 @@ fielddisplay(self,'partitionice','ice partition vector for barystatic contribution'); fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution'); + if isempty(self.external), fielddisplay(self,'external','external solution, of the type solidearthsolution'); end self.settings.disp(); self.surfaceload.disp(); self.lovenumbers.disp(); self.rotational.disp(); + if ~isempty(self.external), + self.external.disp(); + end end % }}} function marshall(self,prefix,md,fid) % {{{ - WriteData(fid,prefix,'object',self,'fieldname','sealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'fieldname','initialsealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','planetradius','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray'); @@ -103,4 +117,7 @@ self.lovenumbers.marshall(prefix,md,fid); self.rotational.marshall(prefix,md,fid); + if ~isempty(self.external), + self.external.marshall(prefix,md,fid); + end %process requested outputs @@ -116,9 +133,12 @@ function savemodeljs(self,fid,modelname) % {{{ - writejs1Darray(fid,[modelname '.solidearth.sealevel'],self.sealevel); + writejs1Darray(fid,[modelname '.solidearth.initialsealevel'],self.initialsealevel); self.settings.savemodeljs(fid,modelname); self.surfaceload.savemodeljs(fid,modelname); self.lovenumbers.savemodeljs(fid,modelname); self.rotational.savemodeljs(fid,modelname); + if ~isempty(self.external), + self.external.savemodeljs(fid,modelname); + end writejscellstring(fid,[modelname '.solidearth.requested_outputs'],self.requested_outputs); writejscellarray(fid,[modelname '.solidearth.transitions'],self.transitions); @@ -126,5 +146,5 @@ end % }}} function self = extrude(self,md) % {{{ - self.sealevel=project3d(md,'vector',self.sealevel,'type','node'); + self.initialsealevel=project3d(md,'vector',self.initialsealevel,'type','node'); end % }}} end Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 25762) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 25763) @@ -14,5 +14,6 @@ ocean_area_scaling = 0; runfrequency = 1; %how many time steps we skip before we run grd_core - computesealevelchange = 0; %will grd_core compute sea level? + computesealevelchange = 1; %will sea-level be coputed? + isgrd = 1; %will GRD patterns be computed? degacc = 0; %degree increment for resolution of Green tables horiz = 0; %compute horizontal deformation @@ -42,5 +43,6 @@ self.rotation=1; self.ocean_area_scaling=0; - self.computesealevelchange=0; + self.isgrd=1; + self.computesealevelchange=1; %numerical discretization accuracy @@ -70,5 +72,5 @@ %a coupler to a planet model is provided. - if self.computesealevelchange, + if self.isgrd, if md.transient.iscoupler, %we are good; @@ -90,5 +92,6 @@ fielddisplay(self,'maxiter','maximum number of nonlinear iterations'); fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); - fielddisplay(self,'computesealevelchange','compute sealevel change from GRD in addition to steric?) default 0'); + fielddisplay(self,'computesealevelchange','compute sealevel change (default 1)'); + fielddisplay(self,'isgrd','compute GRD patterns (default 1)'); fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)'); fielddisplay(self,'rigid','rigid earth graviational potential perturbation'); @@ -110,4 +113,5 @@ WriteData(fid,prefix,'object',self,'fieldname','horiz','name','md.solidearth.settings.horiz','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','computesealevelchange','name','md.solidearth.settings.computesealevelchange','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','isgrd','name','md.solidearth.settings.isgrd','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','glfraction','name','md.solidearth.settings.glfraction','format','Integer'); end % }}} Index: /issm/trunk-jpl/src/m/classes/solidearthsolution.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsolution.m (revision 25762) +++ /issm/trunk-jpl/src/m/classes/solidearthsolution.m (revision 25763) @@ -32,17 +32,12 @@ function md = checkconsistency(self,md,solution,analyses) % {{{ - if ~ismember('SealevelriseAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.isslr==0), - return; - end - md = checkfield(md,'fieldname','solidearth.solution.displacementeast','Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','solidearth.solution.displacementnorth','Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','solidearth.solution.displacementup','Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','solidearth.solution.geoid','Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','solidearth.solution.barystaticsealevel','Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','solidearth.external.displacementeast','Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','solidearth.external.displacementnorth','Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','solidearth.external.displacementup','Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','solidearth.external.geoid','Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','solidearth.external.barystaticsealevel','Inf',1,'timeseries',1); end % }}} function disp(self) % {{{ - disp(sprintf(' solidearth solution:')); - fielddisplay(self,'displacementeast','solid-Earth Eastwards bedrock motion (m/yr)'); fielddisplay(self,'displacementnorth','solid-Earth Northwards bedrock motion (m/yr)'); @@ -53,5 +48,4 @@ end % }}} function marshall(self,prefix,md,fid) % {{{ - WriteData(fid,prefix,'object',self,'data',1,'name','md.solidearth.external.nature','format','Integer'); %code 1 for father class 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); Index: /issm/trunk-jpl/test/NightlyRun/test2002.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 25763) @@ -8,5 +8,5 @@ %solidearth loading: {{{ md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); -md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); +md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1); md.dsl.global_average_thermosteric_sea_level_change=[0;0]; md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); Index: /issm/trunk-jpl/test/NightlyRun/test2002.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2002.py (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2002.py (revision 25763) @@ -19,5 +19,5 @@ #solidearth loading: md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, 1)) -md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, 1)) +md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, 1)) md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1)) md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1)) Index: /issm/trunk-jpl/test/NightlyRun/test2003.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2003.m (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2003.m (revision 25763) @@ -8,5 +8,5 @@ %solidearth loading: {{{ md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); -md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); +md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1); md.dsl.global_average_thermosteric_sea_level_change=[0;0]; md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); Index: /issm/trunk-jpl/test/NightlyRun/test2003.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2003.py (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2003.py (revision 25763) @@ -19,5 +19,5 @@ #solidearth loading: {{{ md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, )) -md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, )) +md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, )) md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, )) md.dsl.sea_surface_height_change_above_geoid=np.zeros((md.mesh.numberofvertices+1, )) Index: /issm/trunk-jpl/test/NightlyRun/test2004.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 25763) @@ -167,5 +167,5 @@ end - md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); + md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1); md.dsl.global_average_thermosteric_sea_level_change=[0;0]; @@ -302,5 +302,5 @@ end - md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); + md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1); md.dsl.global_average_thermosteric_sea_level_change=[0;0]; @@ -349,5 +349,5 @@ sl.transfer('mesh.long'); sl.transfer('solidearth.surfaceload.icethicknesschange'); -sl.transfer('solidearth.sealevel'); +sl.transfer('solidearth.initialsealevel'); sl.transfer('dsl.sea_surface_height_change_above_geoid'); sl.transfer('dsl.sea_water_pressure_change_at_sea_floor'); Index: /issm/trunk-jpl/test/NightlyRun/test2004.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2004.py (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2004.py (revision 25763) @@ -206,5 +206,5 @@ md.solidearth.surfaceload.icethicknesschange = np.mean(delHAIS[md.mesh.elements - 1], axis=1) / 100 - md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, )) + md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, )) md.dsl.global_average_thermostatic_sea_level_change = np.zeros((2, 1)) @@ -335,5 +335,5 @@ md.mask.ocean_levelset[pos] = 1 - md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, )) + md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, )) md.dsl.global_average_thermostatic_sea_level_change = np.zeros((2, 1)) @@ -378,5 +378,5 @@ sl.transfer('mesh.long') sl.transfer('solidearth.surfaceload.icethicknesschange') # -sl.transfer('solidearth.sealevel') +sl.transfer('solidearth.initialsealevel') sl.transfer('dsl.sea_surface_height_change_above_geoid') sl.transfer('dsl.sea_water_pressure_change_at_sea_floor') Index: /issm/trunk-jpl/test/NightlyRun/test2005.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2005.m (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2005.m (revision 25763) @@ -8,5 +8,5 @@ %solidearth loading: {{{ md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); -md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); +md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1); md.dsl.global_average_thermosteric_sea_level_change=[1;0]; md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); Index: /issm/trunk-jpl/test/NightlyRun/test2005.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2005.py (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2005.py (revision 25763) @@ -19,5 +19,5 @@ # Solidearth loading #{{{ md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, 1)) -md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, 1)) +md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, 1)) md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1)) md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1)) Index: /issm/trunk-jpl/test/NightlyRun/test2006.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2006.m (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2006.m (revision 25763) @@ -9,5 +9,5 @@ %solidearth loading: {{{ md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); -md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); +md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1); md.dsl.global_average_thermosteric_sea_level_change=[1;0]; md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); Index: /issm/trunk-jpl/test/NightlyRun/test2006.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2006.py (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2006.py (revision 25763) @@ -22,5 +22,5 @@ # Solidearth loading #{{{ md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, 1)) -md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, 1)) +md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, 1)) md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1)) md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1)) Index: /issm/trunk-jpl/test/NightlyRun/test2010.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2010.m (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2010.m (revision 25763) @@ -15,5 +15,5 @@ md.solidearth.surfaceload.icethicknesschange(pos(6:7))=-1; -md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); +md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1); md.dsl.global_average_thermosteric_sea_level_change=[0;0]; md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); Index: /issm/trunk-jpl/test/NightlyRun/test2010.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2010.py (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2010.py (revision 25763) @@ -27,5 +27,5 @@ md.solidearth.surfaceload.icethicknesschange[pos[5:7]] = -1 -md.solidearth.sealevel = np.zeros(md.mesh.numberofvertices) +md.solidearth.initialsealevel = np.zeros(md.mesh.numberofvertices) md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1)) md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1)) Index: /issm/trunk-jpl/test/NightlyRun/test2424.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2424.m (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2424.m (revision 25763) @@ -20,5 +20,5 @@ md.timestepping.time_step=.1; -md.solidearth.sealevel=newforcing(md.timestepping.start_time, md.timestepping.final_time, md.timestepping.time_step,-200,200,md.mesh.numberofvertices); +md.solidearth.initialsealevel=newforcing(md.timestepping.start_time, md.timestepping.final_time, md.timestepping.time_step,-200,200,md.mesh.numberofvertices); md.cluster=generic('name',oshostname(),'np',3); Index: /issm/trunk-jpl/test/NightlyRun/test2424.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2424.py (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2424.py (revision 25763) @@ -30,5 +30,5 @@ md.timestepping.time_step = .1 -md.solidearth.sealevel = newforcing(md.timestepping.start_time, md.timestepping.final_time, +md.solidearth.initialsealevel = newforcing(md.timestepping.start_time, md.timestepping.final_time, md.timestepping.time_step, -200., 200., md.mesh.numberofvertices) Index: /issm/trunk-jpl/test/NightlyRun/test2425.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2425.m (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2425.m (revision 25763) @@ -31,5 +31,5 @@ md.geometry.bed=md.geometry.bed+1000; md.geometry.surface=md.geometry.surface+1000; -md.solidearth.sealevel=1000*ones(md.mesh.numberofvertices,1); +md.solidearth.initialsealevel=1000*ones(md.mesh.numberofvertices,1); md=solve(md,'Transient','checkconsistency','no'); Index: /issm/trunk-jpl/test/NightlyRun/test2425.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2425.py (revision 25762) +++ /issm/trunk-jpl/test/NightlyRun/test2425.py (revision 25763) @@ -39,5 +39,5 @@ md.geometry.bed = md.geometry.bed + 1000 md.geometry.surface = md.geometry.surface + 1000 -md.solidearth.sealevel = 1000 * np.ones((md.mesh.numberofvertices, )) +md.solidearth.initialsealevel = 1000 * np.ones((md.mesh.numberofvertices, )) md = solve(md, 'Transient', 'checkconsistency', 'no')