Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 19983) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 19984) @@ -111,4 +111,5 @@ ./shared/Numerics/NewtonSolveDnorm.cpp\ ./shared/Numerics/extrema.cpp\ + ./shared/Numerics/legendre.cpp\ ./shared/Numerics/XZvectorsToCoordinateSystem.cpp\ ./shared/Exceptions/Exceptions.cpp\ @@ -464,4 +465,10 @@ endif #}}} +#Slr sources {{{ +if SEALEVELRISE +issm_sources += ./cores/sealevelrise_core.cpp\ + ./analyses/SealevelriseAnalysis.cpp +endif +#}}} #Metis sources {{{ if METIS Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 19983) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 19984) @@ -116,4 +116,7 @@ case LsfReinitializationAnalysisEnum : return new LsfReinitializationAnalysis(); #endif + #ifdef _HAVE_SEALEVELRISE_ + case SealevelriseAnalysisEnum : return new SealevelriseAnalysis(); + #endif default : _error_("enum provided not supported ("<numberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + /*Create inputs: */ + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,MaskOceanLevelsetEnum); + iomodel->FetchDataToInput(elements,SealevelriseDeltathicknessEnum); + +}/*}}}*/ +void SealevelriseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + int nl; + IssmDouble* love_h=NULL; + IssmDouble* love_k=NULL; + + parameters->AddObject(iomodel->CopyConstantObject(SealevelriseReltolEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SealevelriseAbstolEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SealevelriseMaxiterEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SealevelriseRigidEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SealevelriseElasticEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SealevelriseEustaticEnum)); + + iomodel->FetchData(&love_h,&nl,NULL,SealevelriseLoveHEnum); + iomodel->FetchData(&love_k,&nl,NULL,SealevelriseLoveKEnum); + + parameters->AddObject(new DoubleVecParam(SealevelriseLoveHEnum,love_h,nl)); + parameters->AddObject(new DoubleVecParam(SealevelriseLoveKEnum,love_k,nl)); + + /*free ressources: */ + xDelete(love_h); + xDelete(love_k); + +}/*}}}*/ + +/*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_index){/*{{{*/ + _error_("Not implemented yet"); +}/*}}}*/ +void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ + /*Default, do nothing*/ + return; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.h (revision 19984) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.h (revision 19984) @@ -0,0 +1,33 @@ +/*! \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); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,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_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 19983) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 19984) @@ -32,4 +32,5 @@ #include "./MasstransportAnalysis.h" #include "./SmbAnalysis.h" +#include "./SealevelriseAnalysis.h" #include "./MeltingAnalysis.h" #include "./MeshdeformationAnalysis.h" Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 19983) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 19984) @@ -1464,4 +1464,8 @@ bool Element::IsIceInElement(){/*{{{*/ return (this->inputs->Min(MaskIceLevelsetEnum)<0.); +} +/*}}}*/ +bool Element::IsWaterInElement(){/*{{{*/ + return (this->inputs->Max(MaskOceanLevelsetEnum)>0.); } /*}}}*/ @@ -1516,4 +1520,5 @@ name==MaterialsRheologyBbarEnum || name==MaterialsRheologyNEnum || + name==SealevelriseSEnum || name==GiaWEnum || name==GiadWdtEnum || Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 19983) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 19984) @@ -119,4 +119,5 @@ bool IsFloating(); bool IsIceInElement(); + bool IsWaterInElement(); bool IsInput(int name); void LinearFloatingiceMeltingRate(); @@ -299,4 +300,9 @@ virtual void GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y)=0; #endif + #ifdef _HAVE_SEALEVELRISE_ + virtual void Sealevelrise(Vector* pSgi,Vector* pSgo,IssmDouble* peustatic,IssmDouble* Sg_old,IssmDouble* Sgi_old,IssmDouble* Sgo_old,IssmDouble* x,IssmDouble* y,IssmDouble* z,IssmDouble oceanarea)=0; + virtual IssmDouble OceanArea(void)=0; + #endif + }; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 19983) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 19984) @@ -180,4 +180,9 @@ void GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y); #endif + #ifdef _HAVE_SEALEVELRISE_ + void Sealevelrise(Vector* pSgi,Vector* pSgo,IssmDouble* peustatic,IssmDouble* Sg_old,IssmDouble* Sgi_old,IssmDouble* Sgo_old,IssmDouble* x,IssmDouble* y,IssmDouble* z, IssmDouble oceanarea){_error_("not implemented yet!");}; + IssmDouble OceanArea(void){_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 19983) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 19984) @@ -166,4 +166,9 @@ #endif +#ifdef _HAVE_SEALEVELRISE_ + void Sealevelrise(Vector* pSgi,Vector* pSgo,IssmDouble* peustatic,IssmDouble* Sg_old,IssmDouble* Sgi_old,IssmDouble* Sgo_old,IssmDouble* x,IssmDouble* y,IssmDouble* z,IssmDouble oceanarea){_error_("not implemented yet!");}; + IssmDouble OceanArea(void){_error_("not implemented yet!");}; +#endif + #ifdef _HAVE_DAKOTA_ void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 19983) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 19984) @@ -172,4 +172,8 @@ void GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; #endif +#ifdef _HAVE_SEALEVELRISE_ + void Sealevelrise(Vector* pSgi,Vector* pSgo,IssmDouble* peustatic,IssmDouble* Sg_old,IssmDouble* Sgi_old,IssmDouble* Sgo_old,IssmDouble* x,IssmDouble* y,IssmDouble* z,IssmDouble oceanarea){_error_("not implemented yet!");}; + IssmDouble OceanArea(void){_error_("not implemented yet!");}; +#endif #ifdef _HAVE_DAKOTA_ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 19983) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 19984) @@ -1034,4 +1034,23 @@ _assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0); return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2; +} +/*}}}*/ +IssmDouble Tria::GetArea3D(void){/*{{{*/ + + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble x1,y1,z1,x2,y2,z2,x3,y3,z3; + IssmDouble detm1,detm2,detm3; + + /*Get xyz list: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + x1=xyz_list[0][0]; y1=xyz_list[0][1]; z1=xyz_list[0][2]; + x2=xyz_list[1][0]; y2=xyz_list[1][1]; z2=xyz_list[1][2]; + x3=xyz_list[2][0]; y3=xyz_list[2][1]; z3=xyz_list[2][2]; + + detm1=x1*y2 - x2*y1 - x1*y3 + x3*y1 + x2*y3 - x3*y2; + detm2=y1*z2 - y2*z1 - y1*z3 + y3*z1 + y2*z3 - y3*z2; + detm3=x2*z1 - x1*z2 + x1*z3 - x3*z1 - x2*z3 + x3*z2; + + return sqrt(pow(detm1,2) + pow(detm2,2) + pow(detm3,2))/2; } /*}}}*/ @@ -3483,4 +3502,161 @@ #endif +#ifdef _HAVE_SEALEVELRISE_ +void Tria::Sealevelrise(Vector* pSgi,Vector* pSgo,IssmDouble* peustatic,IssmDouble* Sg_old,IssmDouble* Sgi_old,IssmDouble* Sgo_old,IssmDouble* x,IssmDouble* y,IssmDouble* z,IssmDouble oceanarea){ /*{{{*/ + + /*diverse:*/ + int gsize; + bool spherical=true; + IssmDouble x0,y0,z0,a; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble area; + IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) + IssmDouble Me; + IssmDouble rho; + + /*love numbers:*/ + IssmDouble* love_h=NULL; + IssmDouble* love_k=NULL; + IssmDouble* deltalove=NULL; + int nl; + + /*ice properties: */ + IssmDouble rho_ice,rho_water,rho_earth; + + /*other constants: */ + //IssmDouble g; + + /*Solution outputs: */ + Vector* pSolution=NULL; + + /*Initialize eustatic component: do not skip this step :):*/ + IssmDouble eustatic = 0; + + /*Computational flags:*/ + bool computerigid = true; + bool computeelastic= true; + bool computeeustatic = true; + + /*recover material parameters: */ + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); + rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + rho_earth=matpar->GetMaterialParameter(MaterialsEarthDensityEnum); + + /*recover love numbers and computational flags: */ + this->parameters->FindParam(&love_h,&nl,SealevelriseLoveHEnum); + this->parameters->FindParam(&love_k,&nl,SealevelriseLoveKEnum); + this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); + this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); + this->parameters->FindParam(&computeeustatic,SealevelriseEustaticEnum); + + /*how many dofs are we working with here? */ + this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + + /*retrieve coordinates: */ + ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); + + /* Where is the centroid of this element?. To avoid issues with lat,long + * being between [0,180] and [0 360], and issues with jumps at the central + * meridian and poles, we do everything in cartesian coordinate system: */ + x0=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3; + y0=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3; + z0=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3; + + /*radius at this location?:*/ + a=sqrt(pow(xyz_list[0][0],2)+pow(xyz_list[0][1],2)+pow(xyz_list[0][2],2)); //a from Farrel and Clark, Eq 4. + + /*Compute mass of the earth:*/ + Me= rho_earth*4/3*PI*pow(a,3); + + /*Compute area of element:*/ + area=GetArea3D(); + + /*Compute ice thickness or sea level change: */ + if(IsWaterInElement()){ + + IssmDouble phi_I_I_O,phi_O_O_O; + + /*From Sg_old, recover water sea level rise:*/ + I=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; + rho=rho_water; + pSolution=pSgo; + + /*Compute eustatic component: */ + if(computeeustatic){ + phi_I_I_O=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; + phi_O_O_O=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; + eustatic -= (phi_I_I_O + phi_O_O_O); //Watch out the sign "-"! + } + + } + else{ + Input* deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum); + if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); + deltathickness_input->GetInputAverage(&I); + rho=rho_ice; + pSolution=pSgi; + + /*Compute eustatic compoent:*/ + if(computeeustatic)eustatic -= rho_ice*area*I/(oceanarea*rho_water); //Watch out the sign "-"! + } + + + /*Speed up of love number comutation: */ + if(computeelastic){ + deltalove=xNew(nl); + for (int n=0;nSetValue(i,rho*a/Me*I*area*(G_rigid+G_elastic),ADD_VAL); + } + } + + /*Assign output pointer:*/ + *peustatic=eustatic; + + /*Free ressources:*/ + xDelete(love_h); + xDelete(love_k); + xDelete(deltalove); + return; +} +/*}}}*/ +IssmDouble Tria::OceanArea(void){ /*{{{*/ + + if(IsWaterInElement()) return GetArea3D(); + else return 0; + +} +/*}}}*/ +#endif + + #ifdef _HAVE_DAKOTA_ void Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 19983) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 19984) @@ -143,4 +143,8 @@ void GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y); #endif + #ifdef _HAVE_SEALEVELRISE_ + void Sealevelrise(Vector* pSgi,Vector* pSgo,IssmDouble* peustatic,IssmDouble* Sg_old,IssmDouble* Sgi_old,IssmDouble* Sgo_old,IssmDouble* x,IssmDouble* y,IssmDouble* z,IssmDouble oceanarea); + IssmDouble OceanArea(void); + #endif /*}}}*/ /*Tria specific routines:{{{*/ @@ -148,5 +152,6 @@ void AddInput(int input_enum, IssmDouble* values, int interpolation_enum); IssmDouble GetArea(void); - IssmDouble GetAreaIce(void); + IssmDouble GetArea3D(void); + IssmDouble GetAreaIce(void); void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints); void GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 19983) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 19984) @@ -534,4 +534,8 @@ break; + case SealevelriseSolutionEnum: + analyses_temp[numanalyses++]=SealevelriseAnalysisEnum; + break; + case SmbSolutionEnum: analyses_temp[numanalyses++]=SmbAnalysisEnum; @@ -2201,4 +2205,61 @@ /*}}}*/ #endif +#ifdef _HAVE_SEALEVELRISE_ +void FemModel::Sealevelrise(Vector* pSgi,Vector* pSgo, IssmDouble* peustatic, Vector* pSg_old, Vector* pSgi_old,Vector* pSgo_old,IssmDouble* x, IssmDouble* y, IssmDouble* z) { /*{{{*/ + + /*serialized vectors:*/ + IssmDouble* Sg_old=NULL; + IssmDouble* Sgi_old=NULL; + IssmDouble* Sgo_old=NULL; + IssmDouble eustatic=0; + IssmDouble eustatic_cpu=0; + IssmDouble eustatic_cpu_e=0; + IssmDouble oceanarea=0; + IssmDouble oceanarea_cpu=0; + int ns; + + /*Serialize vectors from previous iteration:*/ + + Sg_old=pSg_old->ToMPISerial(); + Sgi_old=pSgi_old->ToMPISerial(); + Sgo_old=pSgo_old->ToMPISerial(); + + /*Go through elements, and add contribution from each element to the deflection vector wg:*/ + ns = elements->Size(); + + /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */ + for(int i=0;i(elements->GetObjectByOffset(i)); + oceanarea_cpu += element->OceanArea(); + } + ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Call the sea level rise core: */ + for(int i=0;i(elements->GetObjectByOffset(i)); + element->Sealevelrise(pSgi,pSgo,&eustatic_cpu_e,Sg_old,Sgi_old,Sgo_old,x,y,z,oceanarea); + eustatic_cpu+=eustatic_cpu_e; + } + if(IssmComm::GetRank()==0)if(VerboseConvergence())_printf_("\n"); + + /*Sum all eustatic components from all cpus:*/ + ISSM_MPI_Reduce (&eustatic_cpu,&eustatic,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&eustatic,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers:*/ + *peustatic=eustatic; + + /*Free ressources:*/ + xDelete(Sg_old); + xDelete(Sgi_old); + xDelete(Sgo_old); +} +/*}}}*/ +#endif + void FemModel::HydrologyEPLupdateDomainx(IssmDouble* pEplcount){ /*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 19983) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 19984) @@ -112,4 +112,7 @@ void Deflection(Vector* wg,Vector* dwgdt, IssmDouble* x, IssmDouble* y); #endif + #ifdef _HAVE_SEALEVELRISE_ + void Sealevelrise(Vector* pSgi,Vector* pSgo, IssmDouble* peustatic, Vector* pSg_old, Vector* pSgi_old,Vector* pSgo_old,IssmDouble* x, IssmDouble* y, IssmDouble* z); + #endif void TimeAdaptx(IssmDouble* pdt); void UpdateConstraintsx(void); Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 19983) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 19984) @@ -58,4 +58,6 @@ mantle_shear_modulus=0; mantle_density=0; + + earth_density=0; poisson=0; @@ -160,4 +162,8 @@ iomodel->Constant(&this->mantle_shear_modulus,MaterialsMantleShearModulusEnum); iomodel->Constant(&this->mantle_density,MaterialsMantleDensityEnum); + + /*slr:*/ + iomodel->Constant(&this->earth_density,MaterialsEarthDensityEnum); + break; default: @@ -255,4 +261,6 @@ matpar->mantle_shear_modulus=this->mantle_shear_modulus; matpar->mantle_density=this->mantle_density; + + matpar->earth_density=this->earth_density; return matpar; @@ -302,4 +310,7 @@ MARSHALLING(mantle_shear_modulus); MARSHALLING(mantle_density); + + //slr: + MARSHALLING(earth_density); //Sea ice: @@ -516,4 +527,5 @@ case MaterialsMantleDensityEnum: return this->mantle_density; case MaterialsMantleShearModulusEnum: return this->mantle_shear_modulus; + case MaterialsEarthDensityEnum: return this->earth_density; default: _error_("Enum "<Data(BaseEnum) && iomodel->Data(ThicknessEnum)); this->sigma = (iomodel->Data(MeshZEnum)[i]-iomodel->Data(BaseEnum)[i])/(iomodel->Data(ThicknessEnum)[i]); + break; + case Domain3DsurfaceEnum: + this->latitute = iomodel->Data(MeshLatEnum)[i]; + this->longitude = iomodel->Data(MeshLongEnum)[i]; + this->R = iomodel->Data(MeshREnum)[i]; break; case Domain2DhorizontalEnum: @@ -121,4 +126,16 @@ IssmDouble Vertex::GetZ(){/*{{{*/ return this->z; +} +/*}}}*/ +IssmDouble Vertex::GetLatitude(){/*{{{*/ + return this->latitute; +} +/*}}}*/ +IssmDouble Vertex::GetLongitude(){/*{{{*/ + return this->longitude; +} +/*}}}*/ +IssmDouble Vertex::GetRadius(){/*{{{*/ + return this->R; } /*}}}*/ @@ -239,11 +256,18 @@ } /*}}}*/ -void Vertex::VertexCoordinates(Vector* vx,Vector* vy,Vector* vz){/*{{{*/ +void Vertex::VertexCoordinates(Vector* vx,Vector* vy,Vector* vz, bool spherical){/*{{{*/ if (this->clone==true) return; - vx->SetValue(this->sid,this->x,INS_VAL); - vy->SetValue(this->sid,this->y,INS_VAL); - vz->SetValue(this->sid,this->z,INS_VAL); + if(!spherical){ + vx->SetValue(this->sid,this->x,INS_VAL); + vy->SetValue(this->sid,this->y,INS_VAL); + vz->SetValue(this->sid,this->z,INS_VAL); + } + else{ + vx->SetValue(this->sid,this->latitute,INS_VAL); + vy->SetValue(this->sid,this->longitude,INS_VAL); + vz->SetValue(this->sid,this->R,INS_VAL); + } return; @@ -252,13 +276,22 @@ /*Methods relating to Vertex, but not internal methods: */ -void GetVerticesCoordinates(IssmDouble* xyz,Vertex** vertices, int numvertices){ /*{{{*/ +void GetVerticesCoordinates(IssmDouble* xyz,Vertex** vertices, int numvertices,bool spherical){ /*{{{*/ _assert_(vertices); _assert_(xyz); - for(int i=0;iGetX(); - xyz[i*3+1]=vertices[i]->GetY(); - xyz[i*3+2]=vertices[i]->GetZ(); + if(!spherical){ + for(int i=0;iGetX(); + xyz[i*3+1]=vertices[i]->GetY(); + xyz[i*3+2]=vertices[i]->GetZ(); + } + } + else{ + for(int i=0;iGetLatitude(); + xyz[i*3+1]=vertices[i]->GetLongitude(); + xyz[i*3+2]=vertices[i]->GetRadius(); + } } }/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Vertex.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Vertex.h (revision 19983) +++ /issm/trunk-jpl/src/c/classes/Vertex.h (revision 19984) @@ -28,4 +28,7 @@ IssmDouble y; IssmDouble z; + IssmDouble latitute; + IssmDouble longitude; + IssmDouble R; IssmDouble sigma; //sigma coordinate: (z-bed)/thickness int connectivity; //number of vertices connected to this vertex @@ -52,4 +55,7 @@ IssmDouble GetY(void); IssmDouble GetZ(void); + IssmDouble GetLatitude(void); + IssmDouble GetLongitude(void); + IssmDouble GetRadius(void); void UpdatePosition(Vector* vx,Vector* vy,Vector* vz,Parameters* parameters,IssmDouble* thickness,IssmDouble* bed); void DistributePids(int* ppidcount); @@ -59,9 +65,9 @@ void SetClone(int* minranks); void ToXYZ(Matrix* matrix); - void VertexCoordinates(Vector* vx,Vector* vy,Vector* vz); + void VertexCoordinates(Vector* vx,Vector* vy,Vector* vz,bool spherical=false); }; /*Methods relating to Vertex object: */ -void GetVerticesCoordinates(IssmDouble* xyz,Vertex** vertices, int numvertices); +void GetVerticesCoordinates(IssmDouble* xyz,Vertex** vertices, int numvertices,bool spherical=false); #endif /* _VERTEX_H */ Index: /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 19983) +++ /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 19984) @@ -59,4 +59,7 @@ solutioncore=&masstransport_core; break; + case SealevelriseSolutionEnum: + solutioncore=&sealevelrise_core; + break; case GiaSolutionEnum: #if _HAVE_GIA_ Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 19983) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 19984) @@ -46,4 +46,5 @@ void dummy_core(FemModel* femmodel); void gia_core(FemModel* femmodel); +void sealevelrise_core(FemModel* femmodel); void smb_core(FemModel* femmodel); void damage_core(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 19984) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 19984) @@ -0,0 +1,173 @@ +/*!\file: sealevelrise_core.cpp + * \brief: core of the SLR solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs); + +void sealevelrise_core(FemModel* femmodel){ + + Vector *Sg = NULL; + Vector *Sg_old = NULL; + Vector *Sgi = NULL; //ice convolution + Vector *Sgi_old = NULL; + Vector *Sgo = NULL; //ocean convolution + Vector *Sgo_old = NULL; + + /*parameters: */ + int count; + bool save_results; + int gsize; + int configuration_type; + bool spherical=true; + bool converged=true; + int max_nonlinear_iterations; + IssmDouble eps_rel; + IssmDouble eps_abs; + IssmDouble *x = NULL; + IssmDouble *y = NULL; + IssmDouble *z = NULL; + IssmDouble eustatic; + + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&max_nonlinear_iterations,SealevelriseMaxiterEnum); + femmodel->parameters->FindParam(&eps_rel,SealevelriseReltolEnum); + femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum); + + if(VerboseSolution()) _printf0_(" computing sea level rise\n"); + + /*Call on core computations: */ + femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum); + + /*first, recover lat,long and radius vectors from vertices: */ + VertexCoordinatesx(&x,&y,&z,femmodel->vertices); //no need for z coordinate + + /*Figure out size of g-set deflection vector and allocate solution vector: */ + gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); + + /*Initialize:*/ + Sg = new Vector(gsize); + Sg->Assemble(); + Sg_old = new Vector(gsize); + Sg_old->Assemble(); + + Sgi = new Vector(gsize); + Sgi->Assemble(); + Sgi_old = new Vector(gsize); + Sgi_old->Assemble(); + + Sgo = new Vector(gsize); + Sgo->Assemble(); + Sgo_old = new Vector(gsize); + Sgo_old->Assemble(); + + count=1; + converged=false; + + /*Start loop: */ + for(;;){ + + //save pointer to old sea level rise + delete Sgi_old; Sgi_old=Sgi; + delete Sgo_old; Sgo_old=Sgo; + delete Sg_old; Sg_old=Sg; + + /*Initialize solution vector: */ + Sg = new Vector(gsize); Sg->Assemble(); + Sgi = new Vector(gsize); Sgi->Assemble(); + Sgo = new Vector(gsize); Sgo->Assemble(); + + /*call the main module: */ + femmodel->Sealevelrise(Sgi,Sgo, &eustatic, Sg_old, Sgi_old, Sgo_old,x, y, z); + + /*assemble solution vectors: */ + Sgi->Assemble(); + Sgo->Assemble(); + + /*Sg is the sum of the ice and ocean convolutions + eustatic component: (eq 4 in Farrel and Clark)*/ + Sgi->Copy(Sg); Sg->AXPY(Sgo,1); + Sg->Shift(eustatic); + + /*convergence criterion:*/ + slrconvergence(&converged,Sg,Sg_old,eps_rel,eps_abs); + + /*Increase count: */ + count++; + if(converged==true){ + break; + } + if(count>=max_nonlinear_iterations){ + _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); + converged=true; + break; + } + + } + if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n"); + + InputUpdateFromVectorx(femmodel,Sg,SealevelriseSEnum,VertexSIdEnum); + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + int outputs[1] = {SealevelriseSEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1); + } + xDelete(x); + xDelete(y); + xDelete(z); + delete Sg_old; + delete Sg; +} + +void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ + + bool converged=true; + IssmDouble ndS,nS; + Vector *dSg = NULL; + + //compute norm(du) and norm(u) if requested + dSg=Sg_old->Duplicate(); Sg_old->Copy(dSg); dSg->AYPX(Sg,-1.0); + ndS=dSg->Norm(NORM_TWO); + + if(!xIsNan(eps_rel)){ + nS=Sg_old->Norm(NORM_TWO); + } + + if (xIsNan(ndS) || xIsNan(nS)) _error_("convergence criterion is NaN!"); + + //clean up + delete dSg; + + //print + if(!xIsNan(eps_rel)){ + if((ndS/nS) " << eps_rel*100 << " %\n"); + converged=false; + } + } + if(!xIsNan(eps_abs)){ + if(ndS " << eps_abs << " \n"); + converged=false; + } + } + + /*assign output*/ + *pconverged=converged; + +} /*}}}*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 19983) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 19984) @@ -95,4 +95,6 @@ /*Fetch data:*/ iomodel->FetchData(6,MeshXEnum,MeshYEnum,MeshZEnum,BaseEnum,ThicknessEnum,MaskIceLevelsetEnum); + if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->FetchData(3,MeshLatEnum,MeshLongEnum,MeshREnum); + CreateNumberNodeToElementConnectivity(iomodel); @@ -103,3 +105,4 @@ /*Free data: */ iomodel->DeleteData(6,MeshXEnum,MeshYEnum,MeshZEnum,BaseEnum,ThicknessEnum,MaskIceLevelsetEnum); + if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->DeleteData(3,MeshLatEnum,MeshLongEnum,MeshREnum); } Index: /issm/trunk-jpl/src/c/modules/VertexCoordinatesx/VertexCoordinatesx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/VertexCoordinatesx/VertexCoordinatesx.cpp (revision 19983) +++ /issm/trunk-jpl/src/c/modules/VertexCoordinatesx/VertexCoordinatesx.cpp (revision 19984) @@ -9,5 +9,5 @@ #include "../../toolkits/toolkits.h" -void VertexCoordinatesx( IssmDouble** px, IssmDouble** py, IssmDouble** pz, Vertices* vertices) { +void VertexCoordinatesx( IssmDouble** px, IssmDouble** py, IssmDouble** pz, Vertices* vertices,bool spherical) { /*output: */ @@ -34,5 +34,5 @@ for(i=0;iSize();i++){ Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i); - vertex->VertexCoordinates(vx,vy,vz); + vertex->VertexCoordinates(vx,vy,vz,spherical); } Index: /issm/trunk-jpl/src/c/modules/VertexCoordinatesx/VertexCoordinatesx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/VertexCoordinatesx/VertexCoordinatesx.h (revision 19983) +++ /issm/trunk-jpl/src/c/modules/VertexCoordinatesx/VertexCoordinatesx.h (revision 19984) @@ -8,5 +8,5 @@ /* local prototypes: */ -void VertexCoordinatesx( IssmDouble** px, IssmDouble** py, IssmDouble** pz, Vertices* vertices); +void VertexCoordinatesx( IssmDouble** px, IssmDouble** py, IssmDouble** pz, Vertices* vertices,bool spherical=false); #endif /* _VERTEX_COORDINATESX_H */ Index: /issm/trunk-jpl/src/m/classes/maskpsl.m =================================================================== --- /issm/trunk-jpl/src/m/classes/maskpsl.m (revision 19983) +++ /issm/trunk-jpl/src/m/classes/maskpsl.m (revision 19984) @@ -64,7 +64,7 @@ end % }}} function marshall(self,md,fid) % {{{ - WriteData(fid,'object',self,'fieldname','groundedice_levelset','format','DoubleMat','mattype',1); - WriteData(fid,'object',self,'fieldname','ice_levelset','format','DoubleMat','mattype',1); - WriteData(fid,'object',self,'fieldname','ocean_levelset','format','DoubleMat','mattype',1); + WriteData(fid,'object',self,'class','mask','fieldname','groundedice_levelset','format','DoubleMat','mattype',1); + WriteData(fid,'object',self,'class','mask','fieldname','ice_levelset','format','DoubleMat','mattype',1); + WriteData(fid,'object',self,'class','mask','fieldname','ocean_levelset','format','DoubleMat','mattype',1); % get mask of vertices of elements with ice Index: /issm/trunk-jpl/src/m/classes/matice.m =================================================================== --- /issm/trunk-jpl/src/m/classes/matice.m (revision 19983) +++ /issm/trunk-jpl/src/m/classes/matice.m (revision 19984) @@ -28,4 +28,7 @@ mantle_density = 0.; + %slr + earth_density = 0; + end methods @@ -100,4 +103,7 @@ self.mantle_density = 3.34; % (g/cm^-3) + %SLR + self.earth_density= 5512; % average density of the Earth, (kg/m^3) + end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -115,4 +121,7 @@ md = checkfield(md,'fieldname','materials.mantle_shear_modulus','>',0,'numel',1); md = checkfield(md,'fieldname','materials.mantle_density','>',0,'numel',1); + end + if ismember(SealevelriseAnalysisEnum(),analyses), + md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1); end @@ -140,4 +149,5 @@ fielddisplay(self,'mantle_shear_modulus','Mantle shear modulus [Pa]'); fielddisplay(self,'mantle_density','Mantle density [g/cm^-3]'); + fielddisplay(self,'earth_density','Mantle density [kg/m^-3]'); end % }}} function marshall(self,md,fid) % {{{ @@ -163,4 +173,5 @@ WriteData(fid,'object',self,'class','materials','fieldname','mantle_shear_modulus','format','Double'); WriteData(fid,'object',self,'class','materials','fieldname','mantle_density','format','Double','scale',10^3); + WriteData(fid,'object',self,'class','materials','fieldname','earth_density','format','Double'); end % }}} function savemodeljs(self,fid,modelname) % {{{ @@ -186,4 +197,5 @@ writejsdouble(fid,[modelname '.materials.mantle_shear_modulus'],self.mantle_shear_modulus); writejsdouble(fid,[modelname '.materials.mantle_density'],self.mantle_density); + writejsdouble(fid,[modelname '.materials.earth_density'],self.earth_density); end % }}} Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 19983) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 19984) @@ -22,4 +22,5 @@ initialization = 0; rifts = 0; + slr = 0; debug = 0; @@ -1064,4 +1065,5 @@ md.friction = friction(); md.rifts = rifts(); + md.slr = slr(); md.timestepping = timestepping(); md.groundingline = groundingline(); @@ -1238,4 +1240,5 @@ disp(sprintf('%19s: %-22s -- %s','initialization' ,['[1x1 ' class(self.initialization) ']'],'initial guess/state')); disp(sprintf('%19s: %-22s -- %s','rifts' ,['[1x1 ' class(self.rifts) ']'],'rifts properties')); + disp(sprintf('%19s: %-22s -- %s','slr' ,['[1x1 ' class(self.slr) ']'],'slr forcings')); disp(sprintf('%19s: %-22s -- %s','debug' ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)')); disp(sprintf('%19s: %-22s -- %s','verbose' ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve')); Index: /issm/trunk-jpl/src/m/classes/slr.m =================================================================== --- /issm/trunk-jpl/src/m/classes/slr.m (revision 19984) +++ /issm/trunk-jpl/src/m/classes/slr.m (revision 19984) @@ -0,0 +1,96 @@ +%SLR class definition +% +% Usage: +% slr=slr(); + +classdef slr + properties (SetAccess=public) + deltathickness = NaN; + maxiter = 0; + reltol = 0; + abstol = 0; + love_h = 0; %provided by PREM model + love_k = 0; %ideam + rigid = 0; + elastic = 0; + eustatic = 0; + end + methods + function self = slr(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + otherwise + error('constructor not supported'); + end + end % }}} + function self = setdefaultparameters(self) % {{{ + + %Convergence criterion: absolute, relative and residual + self.reltol=0.01; %1 percent + self.abstol=0.001; %1 mm of sea level rise + + %maximum of non-linear iterations. + self.maxiter=5; + + %computational flags: + self.rigid=1; + self.elastic=1; + self.eustatic=1; + + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + if ~ismember(SealevelriseAnalysisEnum(),analyses), return; end + md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1); + md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1); + md = checkfield(md,'fieldname','slr.reltol','size',[1 1]); + md = checkfield(md,'fieldname','slr.abstol','size',[1 1]); + md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1); + + %check that love numbers are provided at the same level of accuracy: + if (size(self.love_h,1) ~= size(self.love_k,1)), + error('slr error message: love numbers should be provided at the same level of accuracy'); + end + + end % }}} + function disp(self) % {{{ + disp(sprintf(' slr parameters:')); + + fielddisplay(self,'deltathickness','thickness change (main loading of the slr solution core [m]'); + fielddisplay(self,'reltol','sea level rise relative convergence criterion, NaN: not applied'); + fielddisplay(self,'abstol','sea level rise absolute convergence criterion, NaN: not applied'); + fielddisplay(self,'maxiter','maximum number of nonlinear iterations'); + fielddisplay(self,'love_h','love load number for radial displacement'); + fielddisplay(self,'love_k','love load number for gravitational potential perturbation'); + fielddisplay(self,'rigid','rigid earth graviational potential perturbation'); + fielddisplay(self,'elastic','elastic earth graviational potential perturbation'); + fielddisplay(self,'eustatic','eustatic sea level rise'); + + end % }}} + function marshall(self,md,fid) % {{{ + WriteData(fid,'object',self,'class','sealevelrise','fieldname','deltathickness','format','DoubleMat','mattype',1); + WriteData(fid,'object',self,'class','sealevelrise','fieldname','reltol','format','Double'); + WriteData(fid,'object',self,'class','sealevelrise','fieldname','abstol','format','Double'); + WriteData(fid,'object',self,'class','sealevelrise','fieldname','maxiter','format','Integer'); + WriteData(fid,'object',self,'class','sealevelrise','fieldname','love_h','format','DoubleMat','mattype',1); + WriteData(fid,'object',self,'class','sealevelrise','fieldname','love_k','format','DoubleMat','mattype',1); + WriteData(fid,'object',self,'class','sealevelrise','fieldname','rigid','format','Boolean'); + WriteData(fid,'object',self,'class','sealevelrise','fieldname','elastic','format','Boolean'); + WriteData(fid,'object',self,'class','sealevelrise','fieldname','eustatic','format','Boolean'); + end % }}} + function savemodeljs(self,fid,modelname) % {{{ + + writejs1Darray(fid,[modelname '.srl.deltathickness'],self.deltathickness); + writejsdouble(fid,[modelname '.slr.reltol'],self.reltol); + writejsdouble(fid,[modelname '.slr.abstol'],self.abstol); + writejsdouble(fid,[modelname '.slr.maxiter'],self.maxiter); + writejs1Darray(fid,[modelname '.srl.love_h'],self.love_h); + writejs1Darray(fid,[modelname '.srl.love_k'],self.love_k); + writejsdouble(fid,[modelname '.slr.rigid'],self.rigid); + writejsdouble(fid,[modelname '.slr.eustatic'],self.eustatic); + + end % }}} + end +end Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 19983) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 19984) @@ -85,4 +85,7 @@ analyses=[FlaimAnalysisEnum()]; + case SealevelriseSolutionEnum(), + analyses=[SealevelriseAnalysisEnum()]; + case HydrologySolutionEnum(), analyses=[L2ProjectionBaseAnalysisEnum();HydrologyShreveAnalysisEnum();HydrologyDCInefficientAnalysisEnum();HydrologyDCEfficientAnalysisEnum()];