Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 20998) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 20999) @@ -43,8 +43,13 @@ IssmDouble* love_h=NULL; IssmDouble* love_k=NULL; + IssmDouble* love_l=NULL; bool elastic=false; IssmDouble* G_elastic = NULL; IssmDouble* G_elastic_local = NULL; + IssmDouble* U_elastic = NULL; + IssmDouble* U_elastic_local = NULL; + IssmDouble* H_elastic = NULL; + IssmDouble* H_elastic_local = NULL; int M,m,lower_row,upper_row; IssmDouble degacc=.01; @@ -72,4 +77,5 @@ iomodel->FetchData(&love_h,&nl,NULL,"md.slr.love_h"); iomodel->FetchData(&love_k,&nl,NULL,"md.slr.love_k"); + iomodel->FetchData(&love_l,&nl,NULL,"md.slr.love_l"); /*compute elastic green function for a range of angles*/ @@ -77,4 +83,6 @@ M=reCast(180./degacc+1.); G_elastic=xNew(M); + U_elastic=xNew(M); + H_elastic=xNew(M); /*compute combined legendre + love number (elastic green function:*/ @@ -82,4 +90,6 @@ GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); G_elastic_local=xNew(m); + U_elastic_local=xNew(m); + H_elastic_local=xNew(m); for(int i=lower_row;i(IssmComm::GetSize()); int* displs=xNew(IssmComm::GetSize()); @@ -115,4 +143,6 @@ /*All gather:*/ 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 ressources: */ xDelete(recvcounts); @@ -124,10 +154,19 @@ 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)); /*free ressources: */ xDelete(love_h); xDelete(love_k); + xDelete(love_l); xDelete(G_elastic); xDelete(G_elastic_local); + xDelete(U_elastic); + xDelete(U_elastic_local); + xDelete(H_elastic); + xDelete(H_elastic_local); } Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 20998) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 20999) @@ -1501,4 +1501,8 @@ name==MaterialsRheologyEsbarEnum || name==SealevelEnum || + name==SealevelUmotionEnum || + name==SealevelNmotionEnum || + name==SealevelEmotionEnum || + name==SealevelAbsoluteEnum || name==SealevelEustaticEnum || name==SealevelriseDeltathicknessEnum || Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 20998) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 20999) @@ -303,4 +303,5 @@ virtual void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0; virtual void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0; + virtual void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 20998) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 20999) @@ -188,4 +188,5 @@ void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea){_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 20998) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 20999) @@ -172,4 +172,5 @@ void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; IssmDouble OceanArea(void){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 20998) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 20999) @@ -178,4 +178,5 @@ void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; IssmDouble OceanArea(void){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20998) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20999) @@ -3883,6 +3883,200 @@ } /*}}}*/ +void Tria::SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ + + /*diverse:*/ + int gsize; + bool spherical=true; + IssmDouble llr_list[NUMVERTICES][3]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble area; + IssmDouble I; //change in relative sea level or ice thickness + IssmDouble late,longe,re; + IssmDouble lati,longi,ri; + IssmDouble rho_ice,rho_water,rho_earth; + IssmDouble minlong=400; + IssmDouble maxlong=-20; + + /*precomputed elastic green functions:*/ + IssmDouble* U_elastic_precomputed = NULL; + IssmDouble* H_elastic_precomputed = NULL; + int M; + + /*computation of Green functions:*/ + IssmDouble* U_elastic= NULL; + IssmDouble* N_elastic= NULL; + IssmDouble* E_elastic= NULL; + + /*optimization:*/ + bool store_green_functions=false; + + /*computational flags:*/ + bool computerigid = true; + bool computeelastic= true; + + /*early return if we are not on the ocean or on an ice cap:*/ + if(!(this->inputs->Max(MaskIceLevelsetEnum)<0) && !IsWaterInElement()) return; + + /*recover computational flags: */ + this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); + this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); + + /*early return if rigid or elastic not requested:*/ + if(!computerigid && !computeelastic) return; + + /*recover material parameters: */ + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); + rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + rho_earth=matpar->GetMaterialParameter(MaterialsEarthDensityEnum); + + /*how many dofs are we working with here? */ + this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + + /*compute area of element:*/ + area=GetAreaSpherical(); + + /*element centroid (spherical): */ + /* Where is the centroid of this element?:{{{*/ + ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical); + + minlong=400; maxlong=-20; + for (int i=0;imaxlong)maxlong=llr_list[i][1]; + if(llr_list[i][1]180){ + if (llr_list[0][1]==0)llr_list[0][1]=360; + if (llr_list[1][1]==0)llr_list[1][1]=360; + if (llr_list[2][1]==0)llr_list[2][1]=360; + } + + // correction at the north pole + if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; + if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0; + if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0; + + //correction at the south pole + if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; + if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0; + if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0; + + late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0; + longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0; + + late=90-late; + if(longe>180)longe=(longe-180)-180; + + late=late/180*PI; + longe=longe/180*PI; + /*}}}*/ + + /*figure out gravity center of our element (Cartesian): */ + IssmDouble x_element, y_element, z_element; + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; + y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; + z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0; + + if(computeelastic){ + + /*recover elastic Green's functions for displacement:*/ + DoubleVecParam* U_parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter); + DoubleVecParam* H_parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter); + U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M); + H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M); + + /*initialize: */ + U_elastic=xNewZeroInit(gsize); + N_elastic=xNewZeroInit(gsize); + E_elastic=xNewZeroInit(gsize); + } + + int* indices=xNew(gsize); + IssmDouble* U_values=xNewZeroInit(gsize); + IssmDouble* N_values=xNewZeroInit(gsize); + IssmDouble* E_values=xNewZeroInit(gsize); + + for(int i=0;i(alpha/PI*(M-1)); + U_elastic[i] += U_elastic_precomputed[index]; + N_elastic[i] += H_elastic_precomputed[index]*N_azim; + E_elastic[i] += H_elastic_precomputed[index]*E_azim; + } + + /*Add all components to the pUp solution vectors:*/ + if(computerigid){ + U_values[i]+=0; N_values[i]+=0; E_values[i]+=0; + } + if(computeelastic){ + + if(this->inputs->Max(MaskIceLevelsetEnum)<0){ + + /*Compute ice thickness change: */ + Input* deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum); + if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); + deltathickness_input->GetInputAverage(&I); + + U_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i]; + N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i]; + E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i]; + } + else if(IsWaterInElement()) { + + /*From Sg, recover water sea level rise:*/ + I=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; + + U_values[i]+=3*rho_water/rho_earth*area/eartharea*I*U_elastic[i]; + N_values[i]+=3*rho_water/rho_earth*area/eartharea*I*N_elastic[i]; + E_values[i]+=3*rho_water/rho_earth*area/eartharea*I*E_elastic[i]; + } + } + } + pUp->SetValues(gsize,indices,U_values,ADD_VAL); + pNorth->SetValues(gsize,indices,N_values,ADD_VAL); + pEast->SetValues(gsize,indices,E_values,ADD_VAL); + + /*free ressources:*/ + xDelete(indices); + xDelete(U_values); xDelete(N_values); xDelete(E_values); + + /*Free ressources:*/ + if(computeelastic) { + xDelete(U_elastic); xDelete(N_elastic); xDelete(E_elastic); + } + + return; +} +/*}}}*/ #endif - #ifdef _HAVE_DAKOTA_ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 20998) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 20999) @@ -149,4 +149,5 @@ void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); + void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea); #endif /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 20998) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 20999) @@ -2399,4 +2399,66 @@ } /*}}}*/ +void FemModel::SealevelriseGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/ + + /*serialized vectors:*/ + IssmDouble* Sg=NULL; + + IssmDouble oceanarea=0; + IssmDouble oceanarea_cpu=0; + IssmDouble eartharea=0; + IssmDouble eartharea_cpu=0; + + int ns,nsmax; + + /*Serialize vectors from previous iteration:*/ + Sg=pSg->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(); + eartharea_cpu += element->GetAreaSpherical(); + } + 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()); + + ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Figure out max of ns: */ + ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + + /*Call the sea level rise core: */ + for(int i=0;i(elements->GetObjectByOffset(i)); + element->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,oceanarea,eartharea); + } + if(i%100==0){ + pUp->Assemble(); + pNorth->Assemble(); + pEast->Assemble(); + } + } + + /*One last time: */ + pUp->Assemble(); + pNorth->Assemble(); + pEast->Assemble(); + + /*Free ressources:*/ + xDelete(Sg); + xDelete(latitude); + xDelete(longitude); + xDelete(radius); + xDelete(xx); + xDelete(yy); + xDelete(zz); +} +/*}}}*/ IssmDouble FemModel::SealevelriseOceanAverage(Vector* Sg) { /*{{{*/ @@ -2423,5 +2485,4 @@ ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - /*Free ressources:*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 20998) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 20999) @@ -116,4 +116,5 @@ void SealevelriseEustatic(Vector* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution); + void SealevelriseGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); IssmDouble SealevelriseOceanAverage(Vector* Sg); #endif Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 20998) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 20999) @@ -13,5 +13,9 @@ Vector *Sg = NULL; - Vector *Sg_eustatic = NULL; + Vector *Sg_absolute = NULL; + Vector *Sg_eustatic = NULL; + Vector *U_radial = NULL; + Vector *U_north = NULL; + Vector *U_east = NULL; bool save_results,isslr,iscoupler; int configuration_type; @@ -20,4 +24,14 @@ char **requested_outputs = NULL; + /*additional parameters: */ + int gsize; + bool spherical=true; + IssmDouble *latitude = NULL; + IssmDouble *longitude = NULL; + IssmDouble *radius = NULL; + IssmDouble *xx = NULL; + IssmDouble *yy = NULL; + IssmDouble *zz = NULL; + /*Recover some parameters: */ femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); @@ -26,4 +40,13 @@ femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); + + /*first, recover lat,long and radius vectors from vertices: */ + VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); + + /*recover x,y,z vectors from vertices: */ + VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); + + /*Figure out size of g-set deflection vector and allocate solution vector: */ + gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); /*several cases here, depending on value of iscoupler and isslr: @@ -57,8 +80,28 @@ Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //ocean loading tems (2nd and 5th terms on the RHS of Farrel and Clark) - + /*get results into elements:*/ - InputUpdateFromSolutionx(femmodel,Sg); - + //InputUpdateFromSolutionx(femmodel,Sg); // from Eric + InputUpdateFromVectorx(femmodel,Sg,SealevelEnum,VertexSIdEnum); + + /*compute other geodetic signatures, such as absolute sea level chagne, components of 3-D crustal motion: */ + /*Initialize:*/ + U_radial = new Vector(gsize); + U_north = new Vector(gsize); + U_east = new Vector(gsize); + Sg_absolute = new Vector(gsize); + + /*call the geodetic main modlule:*/ + femmodel->SealevelriseGeodetic(U_radial,U_north,U_east,Sg,latitude,longitude,radius,xx,yy,zz); + + /*compute: absolute sea level change = relative sea level change + vertical motion*/ + Sg->Copy(Sg_absolute); Sg_absolute->AXPY(U_radial,1); + + /*get results into elements:*/ + InputUpdateFromVectorx(femmodel,U_radial,SealevelUmotionEnum,VertexSIdEnum); // radial displacement + InputUpdateFromVectorx(femmodel,U_north,SealevelNmotionEnum,VertexSIdEnum); // north motion + InputUpdateFromVectorx(femmodel,U_east,SealevelEmotionEnum,VertexSIdEnum); // east motion + InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelAbsoluteEnum,VertexSIdEnum); + if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); @@ -66,5 +109,5 @@ femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); } - + if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx(); @@ -72,4 +115,8 @@ delete Sg; delete Sg_eustatic; + delete U_radial; + delete U_north; + delete U_east; + delete Sg_absolute; if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} } Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core_eustatic.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core_eustatic.cpp (revision 20998) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core_eustatic.cpp (revision 20999) @@ -58,2 +58,3 @@ return Sgi; } + Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp (revision 20998) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp (revision 20999) @@ -12,5 +12,5 @@ void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs); -Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* Sg_eustatic){ +Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* Sg_eustatic){ /*{{{*/ Vector *Sg = NULL; @@ -35,5 +35,4 @@ IssmDouble *radius = NULL; IssmDouble eustatic; - /*Recover some parameters: */ @@ -112,5 +111,5 @@ return Sg; -} +} /*}}}*/ void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ @@ -157,2 +156,3 @@ } /*}}}*/ + Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 20998) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 20999) @@ -1065,4 +1065,8 @@ SealevelriseAnalysisEnum, SealevelEnum, + SealevelUmotionEnum, + SealevelNmotionEnum, + SealevelEmotionEnum, + SealevelAbsoluteEnum, SealevelEustaticEnum, SealevelriseDeltathicknessEnum, @@ -1078,4 +1082,6 @@ SealevelriseRotationEnum, SealevelriseGElasticEnum, + SealevelriseUElasticEnum, + SealevelriseHElasticEnum, SealevelriseDegaccEnum, SealevelriseTransitionsEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 20998) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 20999) @@ -1017,4 +1017,8 @@ case SealevelriseAnalysisEnum : return "SealevelriseAnalysis"; case SealevelEnum : return "Sealevel"; + case SealevelUmotionEnum : return "SealevelUmotion"; + case SealevelNmotionEnum : return "SealevelNmotion"; + case SealevelEmotionEnum : return "SealevelEmotion"; + case SealevelAbsoluteEnum : return "SealevelAbsolute"; case SealevelEustaticEnum : return "SealevelEustatic"; case SealevelriseDeltathicknessEnum : return "SealevelriseDeltathickness"; @@ -1030,4 +1034,6 @@ case SealevelriseRotationEnum : return "SealevelriseRotation"; case SealevelriseGElasticEnum : return "SealevelriseGElastic"; + case SealevelriseUElasticEnum : return "SealevelriseUElastic"; + case SealevelriseHElasticEnum : return "SealevelriseHElastic"; case SealevelriseDegaccEnum : return "SealevelriseDegacc"; case SealevelriseTransitionsEnum : return "SealevelriseTransitions"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 20998) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 20999) @@ -1041,4 +1041,8 @@ else if (strcmp(name,"SealevelriseAnalysis")==0) return SealevelriseAnalysisEnum; else if (strcmp(name,"Sealevel")==0) return SealevelEnum; + else if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum; + else if (strcmp(name,"SealevelNmotion")==0) return SealevelNmotionEnum; + else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum; + else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum; else if (strcmp(name,"SealevelEustatic")==0) return SealevelEustaticEnum; else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum; @@ -1054,4 +1058,6 @@ else if (strcmp(name,"SealevelriseRotation")==0) return SealevelriseRotationEnum; else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum; + else if (strcmp(name,"SealevelriseUElastic")==0) return SealevelriseUElasticEnum; + else if (strcmp(name,"SealevelriseHElastic")==0) return SealevelriseHElasticEnum; else if (strcmp(name,"SealevelriseDegacc")==0) return SealevelriseDegaccEnum; else if (strcmp(name,"SealevelriseTransitions")==0) return SealevelriseTransitionsEnum;