Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 20006) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 20007) @@ -468,5 +468,7 @@ if SEALEVELRISE issm_sources += ./cores/sealevelrise_core.cpp\ - ./analyses/SealevelriseAnalysis.cpp + ./cores/sealevelrise_core_eustatic.cpp\ + ./cores/sealevelrise_core_noneustatic.cpp\ + ./analyses/SealevelriseAnalysis.cpp endif #}}} Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 20006) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 20007) @@ -302,6 +302,9 @@ #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 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 IssmDouble OceanAverage(IssmDouble* Sg)=0; virtual IssmDouble OceanArea(void)=0; + virtual IssmDouble GetArea3D(void)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 20006) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 20007) @@ -171,4 +171,5 @@ void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); + IssmDouble GetArea3D(void){_error_("not implemented yet!");}; #ifdef _HAVE_DAKOTA_ @@ -181,6 +182,8 @@ #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!");}; + 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!");}; IssmDouble OceanArea(void){_error_("not implemented yet!");}; + IssmDouble OceanAverage(IssmDouble* Sg){_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 20006) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 20007) @@ -161,4 +161,5 @@ void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");}; void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");}; + IssmDouble GetArea3D(void){_error_("not implemented yet!");}; #ifdef _HAVE_GIA_ @@ -167,6 +168,8 @@ #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!");}; + 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!");}; IssmDouble OceanArea(void){_error_("not implemented yet!");}; + IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 20006) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 20007) @@ -168,4 +168,5 @@ void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); + IssmDouble GetArea3D(void){_error_("not implemented yet!");}; #ifdef _HAVE_GIA_ @@ -173,6 +174,8 @@ #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!");}; + 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!");}; IssmDouble OceanArea(void){_error_("not implemented yet!");}; + IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20006) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20007) @@ -3503,16 +3503,16 @@ #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){ /*{{{*/ +void Tria::SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ /*diverse:*/ int gsize; bool spherical=true; - IssmDouble x0,y0,z0,a; - IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble llr_list[NUMVERTICES][3]; IssmDouble area; IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) - IssmDouble Me; IssmDouble rho; - + IssmDouble late,longe,re; + IssmDouble lati,longi,ri; + /*love numbers:*/ IssmDouble* love_h=NULL; @@ -3529,7 +3529,4 @@ IssmDouble rho_ice,rho_water,rho_earth; - /*Solution outputs: */ - Vector* pSolution=NULL; - /*Initialize eustatic component: do not skip this step :):*/ IssmDouble eustatic = 0; @@ -3539,5 +3536,12 @@ bool computeelastic= true; bool computeeustatic = true; + + /*early return if we are not on an ice cap:*/ + if(IsWaterInElement() | !IsIceInElement()){ + *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! + return; + } + /*recover material parameters: */ rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); @@ -3554,5 +3558,5 @@ /*recover legendre coefficients if they have been precomputed:*/ this->parameters->FindParam(&legendre_precompute,SealevelriseLegendrePrecomputeEnum); - this->parameters->FindParam(&legendre_coefficients,&M,NULL,SealevelriseLegendreCoefficientsEnum); + if(legendre_precompute)this->parameters->FindParam(&legendre_coefficients,&M,NULL,SealevelriseLegendreCoefficientsEnum); /*how many dofs are we working with here? */ @@ -3560,54 +3564,50 @@ /*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); + ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical); + + /* Where is the centroid of this element?:*/ + IssmDouble minlong=400; + IssmDouble 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; /*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 if(IsIceInElement()){ - 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 "-"! - } - else if(IsLandInElement()){ - //do nothing - } - else _error_("Tria::Sealevelrise error: partition of land, ice and ocean is not complete!"); + /*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); + + /*Compute eustatic compoent:*/ + if(computeeustatic) eustatic += rho_ice*area*I/(oceanarea*rho_water); /*Speed up of love number comutation: */ @@ -3625,11 +3625,11 @@ IssmDouble Pn1,Pn2,Pn; //first two legendre polynomials IssmDouble cosalpha; - - /*Compute alpha angle between centroid and current vertex and cosine of this angle: */ - alpha = acos( - (x[i]*x0+y[i]*y0+z[i]*z0) //scalar product of two position vectors - / sqrt(pow(x0,2.0)+pow(y0,2.0)+pow(z0,2.0)) //norm of first position vector - / sqrt(pow(x[i],2.0)+pow(y[i],2.0)+pow(z[i],2.0)) //norm of second position vector - ); + IssmDouble delPhi,delLambda; + + /*Compute alpha angle between centroid and current vertex : */ + lati=latitude[i]/180*PI; longi=longitude[i]/180*PI; + + delPhi=fabs(lati-late); delLambda=fabs(longi-longe); + alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); cosalpha=cos(alpha); //compute this to speed up the elastic component. @@ -3652,5 +3652,6 @@ /*Add all components to the pSgi or pSgo solution vectors:*/ - pSolution->SetValue(i,rho*a/Me*I*area*(G_rigid+G_elastic),ADD_VAL); + pSgi->SetValue(i,3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic),ADD_VAL); + } } @@ -3660,8 +3661,219 @@ /*Free ressources:*/ - xDelete(love_h); - xDelete(love_k); - xDelete(deltalove); + if(computeelastic){ + xDelete(love_h); + xDelete(love_k); + xDelete(deltalove); + } return; +} +/*}}}*/ +void Tria::SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ + + /*diverse:*/ + int gsize; + bool spherical=true; + IssmDouble llr_list[NUMVERTICES][3]; + IssmDouble area; + IssmDouble I; //change in water water level(Farrel and Clarke, Equ. 4) + IssmDouble late,longe,re; + IssmDouble lati,longi,ri; + IssmDouble minlong=400; + IssmDouble maxlong=-20; + + /*love numbers:*/ + IssmDouble* love_h=NULL; + IssmDouble* love_k=NULL; + IssmDouble* deltalove=NULL; + int nl; + + /*legendre coefficients:*/ + bool legendre_precompute=false; + IssmDouble* legendre_coefficients=NULL; + int M; + + /*ice properties: */ + IssmDouble rho_ice,rho_water,rho_earth; + + /*Computational flags:*/ + bool computerigid = true; + bool computeelastic= true; + bool computeeustatic = true; + + /*G_rigid and G_elasti variables, to speed up the computations: */ + DoubleArrayInput* G_rigid_input=NULL; + IssmDouble* G_rigid = NULL; int G_rigid_M; + bool compute_G_rigid=false; + DoubleArrayInput* G_elastic_input=NULL; + IssmDouble* G_elastic = NULL; int G_elastic_M; + bool compute_G_elastic=false; + + /*early return if we are not on the ocean:*/ + if (!IsWaterInElement())return; + + /*recover computational flags: */ + this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); + this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); + this->parameters->FindParam(&computeeustatic,SealevelriseEustaticEnum); + + /*early return if rigid or elastic not requested:*/ + if(!computerigid && !computeelastic) return; + + /*Try and recover, if it exists, G_rigid and G_elastic:*/ + G_rigid_input=(DoubleArrayInput*)this->inputs->GetInput(SealevelriseGRigidEnum); + if(G_rigid_input){ + compute_G_rigid=false; + G_rigid_input->GetValues(&G_rigid,&G_rigid_M); + } + else if(computerigid)compute_G_rigid=true; + + G_elastic_input=(DoubleArrayInput*)this->inputs->GetInput(SealevelriseGElasticEnum); + if(G_elastic_input){ + compute_G_elastic=false; + G_elastic_input->GetValues(&G_elastic,&G_elastic_M); + } + else if(computeelastic)compute_G_elastic=true; + + /*recover material parameters: */ + rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + rho_earth=matpar->GetMaterialParameter(MaterialsEarthDensityEnum); + + /*how many dofs are we working with here? */ + this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + + /*From Sg_old, recover water sea level rise:*/ + I=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; + + /*Compute area of element:*/ + area=GetArea3D(); + + /* 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; + /*}}}*/ + if(compute_G_elastic){ + /*recover love numbers and legendre_coefficients:{{{*/ + this->parameters->FindParam(&love_h,&nl,SealevelriseLoveHEnum); + this->parameters->FindParam(&love_k,&nl,SealevelriseLoveKEnum); + + /*recover legendre coefficients if they have been precomputed:*/ + this->parameters->FindParam(&legendre_precompute,SealevelriseLegendrePrecomputeEnum); + if(legendre_precompute)this->parameters->FindParam(&legendre_coefficients,&M,NULL,SealevelriseLegendreCoefficientsEnum); + + /*Speed up of love number comutation for elastic mode: */ + deltalove=xNew(nl); + for (int n=0;n(gsize); + } + if(compute_G_rigid){ + /*Initialize G_rigid and G_elastic:*/ + G_rigid=xNew(gsize); + } + + for(int i=0;i((M-1)*(cosalpha+1.0)/2.0)*nl+n]; + } + else{ + for(int n=0;n(cosalpha),n); Pn1=Pn2; Pn2=Pn; + G_elastic[i] += deltalove[n]*Pn; + } + } + } + } + + /*Add all components to the pSgi or pSgo solution vectors:*/ + if(computerigid)pSgo->SetValue(i,3*rho_water/rho_earth*area/eartharea*I*G_rigid[i],ADD_VAL); + if(computeelastic)pSgo->SetValue(i,3*rho_water/rho_earth*area/eartharea*I*G_elastic[i],ADD_VAL); + } + + /*Save G_rigid and G_elastic if computed:*/ + if(G_rigid)this->inputs->AddInput(new DoubleArrayInput(SealevelriseGRigidEnum,G_rigid,gsize)); + if(G_elastic)this->inputs->AddInput(new DoubleArrayInput(SealevelriseGElasticEnum,G_elastic,gsize)); + + /*Free ressources:*/ + if(compute_G_elastic){ + xDelete(love_h); + xDelete(love_k); + xDelete(deltalove); + xDelete(G_elastic); + } + if(compute_G_rigid) xDelete(G_rigid); + + return; +} +/*}}}*/ +IssmDouble Tria::OceanAverage(IssmDouble* Sg){ /*{{{*/ + + if(IsWaterInElement()){ + + IssmDouble area; + + /*Compute area of element:*/ + area=GetArea3D(); + + /*Average Sg over vertices:*/ + IssmDouble Sg_avg=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; + + /*return: */ + return area*Sg_avg; + } + else return 0; + } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 20006) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 20007) @@ -144,5 +144,7 @@ #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); + 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); + IssmDouble OceanAverage(IssmDouble* Sg); IssmDouble OceanArea(void); #endif Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 20006) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 20007) @@ -2206,10 +2206,7 @@ #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) { /*{{{*/ +void FemModel::SealevelriseEustatic(Vector* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/ /*serialized vectors:*/ - IssmDouble* Sg_old=NULL; - IssmDouble* Sgi_old=NULL; - IssmDouble* Sgo_old=NULL; IssmDouble eustatic=0; IssmDouble eustatic_cpu=0; @@ -2217,12 +2214,8 @@ IssmDouble oceanarea=0; IssmDouble oceanarea_cpu=0; + IssmDouble eartharea=0; + IssmDouble eartharea_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(); @@ -2232,7 +2225,11 @@ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); oceanarea_cpu += element->OceanArea(); + eartharea_cpu+= element->GetArea3D(); } 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()); /*Call the sea level rise core: */ @@ -2242,5 +2239,5 @@ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->Sealevelrise(pSgi,pSgo,&eustatic_cpu_e,Sg_old,Sgi_old,Sgo_old,x,y,z,oceanarea); + element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea); eustatic_cpu+=eustatic_cpu_e; } @@ -2254,8 +2251,81 @@ *peustatic=eustatic; +} +/*}}}*/ +void FemModel::SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution){/*{{{*/ + + /*serialized vectors:*/ + IssmDouble* Sg_old=NULL; + + IssmDouble oceanarea=0; + IssmDouble oceanarea_cpu=0; + IssmDouble eartharea=0; + IssmDouble eartharea_cpu=0; + + int ns; + + /*Serialize vectors from previous iteration:*/ + Sg_old=pSg_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(); + eartharea_cpu+= element->GetArea3D(); + } + 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()); + + /*Call the sea level rise core: */ + for(int i=0;i(elements->GetObjectByOffset(i)); + element->SealevelriseNonEustatic(pSgo, Sg_old,latitude,longitude,radius,oceanarea,eartharea); + } + if(verboseconvolution)if(IssmComm::GetRank()==0)if(VerboseConvergence())_printf_("\n"); + + /*Free ressources:*/ xDelete(Sg_old); - xDelete(Sgi_old); - xDelete(Sgo_old); +} +/*}}}*/ +IssmDouble FemModel::SealevelriseOceanAverage(Vector* Sg) { /*{{{*/ + + IssmDouble* Sg_serial=NULL; + IssmDouble oceanvalue,oceanvalue_cpu; + IssmDouble oceanarea,oceanarea_cpu; + + /*Serialize vectors from previous iteration:*/ + Sg_serial=Sg->ToMPISerial(); + + /*Initialize:*/ + oceanvalue_cpu=0; + oceanarea_cpu=0; + + /*Go through elements, and add contribution from each element and divide by overall ocean area:*/ + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + oceanarea_cpu += element->OceanArea(); + oceanvalue_cpu += element->OceanAverage(Sg_serial); + } + 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 (&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:*/ + xDelete(Sg_serial); + + return oceanvalue/oceanarea; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 20006) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 20007) @@ -113,5 +113,7 @@ #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); + 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); + IssmDouble SealevelriseOceanAverage(Vector* Sg); #endif void TimeAdaptx(IssmDouble* pdt); Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 20006) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 20007) @@ -46,7 +46,9 @@ 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); +void sealevelrise_core(FemModel* femmodel); +Vector * sealevelrise_core_eustatic(FemModel* femmodel); +Vector * sealevelrise_core_noneustatic(FemModel* femmodel,Vector* Sg_eustatic); IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel); Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 20006) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 20007) @@ -10,109 +10,26 @@ #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; + Vector *Sg_eustatic = NULL; 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); + int configuration_type; if(VerboseSolution()) _printf0_(" computing sea level rise\n"); - /*Call on core computations: */ + /*Recover some parameters: */ + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + + /*set configuration: */ femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum); - /*first, recover lat,long and radius vectors from vertices: */ - VertexCoordinatesx(&x,&y,&z,femmodel->vertices); //no need for z coordinate + /*call two sub cores:*/ + Sg_eustatic=sealevelrise_core_eustatic(femmodel); //generalized eustatic (Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS. - /*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(); + Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //sea-level rise dependent terms (2nd and 5th terms on the RHS) - 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"); - + /*get results out:*/ InputUpdateFromVectorx(femmodel,Sg,SealevelriseSEnum,VertexSIdEnum); @@ -122,52 +39,5 @@ femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1); } - xDelete(x); - xDelete(y); - xDelete(z); - delete Sg_old; delete Sg; + delete Sg_eustatic; } - -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/cores/sealevelrise_core_eustatic.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core_eustatic.cpp (revision 20007) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core_eustatic.cpp (revision 20007) @@ -0,0 +1,60 @@ +/*!\file: sealevelrise_core_eustatic.cpp + * \brief: eustatic core of the SLR solution (terms that are constant with respect to sea-level) + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +Vector* sealevelrise_core_eustatic(FemModel* femmodel){ + + Vector *Sgi = NULL; + IssmDouble Sgi_oceanaverage = 0; + + /*parameters: */ + int configuration_type; + int gsize; + bool spherical=true; + IssmDouble *latitude = NULL; + IssmDouble *longitude = NULL; + IssmDouble *radius = NULL; + + /*outputs:*/ + IssmDouble eustatic=0; + + /*recover parameters:*/ + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + + /*first, recover lat,long and radius vectors from vertices: */ + VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); + + /*Figure out size of g-set deflection vector and allocate solution vector: */ + gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); + + /*Initialize:*/ + Sgi = new Vector(gsize); + Sgi->Assemble(); + + /*call the eustatic main module: */ + femmodel->SealevelriseEustatic(Sgi,&eustatic, latitude, longitude, radius); //this computes + + /*assemble solution vector: */ + Sgi->Assemble(); + + /*we need to average Sgi over the ocean: RHS term 4 in Eq.4 of Farrel and clarke. Only the elements can do that: */ + Sgi_oceanaverage=femmodel->SealevelriseOceanAverage(Sgi); + + /*Sg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the + * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/ + Sgi->Shift(-eustatic-Sgi_oceanaverage); + + xDelete(latitude); + xDelete(longitude); + xDelete(radius); + + /*return:*/ + return Sgi; +} Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp (revision 20007) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp (revision 20007) @@ -0,0 +1,158 @@ +/*!\file: sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5. + * \brief: non eustatic core of the SLR solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs); + +Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* Sg_eustatic){ + + Vector *Sg = NULL; + Vector *Sg_old = NULL; + + Vector *Sgo = NULL; //ocean convolution of the perturbation to gravity potential. + IssmDouble Sgo_oceanaverage = 0; //average of Sgo over the ocean. + + /*parameters: */ + int count; + bool save_results; + int gsize; + int configuration_type; + bool spherical=true; + bool converged=true; + bool verboseconvolution=true; + int max_nonlinear_iterations; + IssmDouble eps_rel; + IssmDouble eps_abs; + IssmDouble *latitude = NULL; + IssmDouble *longitude = NULL; + IssmDouble *radius = NULL; + IssmDouble eustatic; + + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&max_nonlinear_iterations,SealevelriseMaxiterEnum); + femmodel->parameters->FindParam(&eps_rel,SealevelriseReltolEnum); + femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum); + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + + /*first, recover lat,long and radius vectors from vertices: */ + VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); + + /*Figure out size of g-set deflection vector and allocate solution vector: */ + gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); + + /*Initialize:*/ + Sg = new Vector(gsize); + Sg->Assemble(); + Sg_eustatic->Copy(Sg); //first initialize Sg with the eustatic component computed in sealevelrise_core_eustatic. + + Sg_old = new Vector(gsize); + Sg_old->Assemble(); + + Sgo = new Vector(gsize); + Sgo->Assemble(); + + count=1; + converged=false; + + /*Start loop: */ + for(;;){ + + //save pointer to old sea level rise + delete Sg_old; Sg_old=Sg; + + /*Initialize solution vector: */ + Sg = new Vector(gsize); Sg->Assemble(); + Sgo = new Vector(gsize); Sgo->Assemble(); + + /*call the non eustatic module: */ + femmodel->SealevelriseNonEustatic(Sgo, Sg_old, latitude, longitude, radius,verboseconvolution); + + /*assemble solution vector: */ + Sgo->Assemble(); + + /*we need to average Sgo over the ocean: RHS term 5 in Eq.4 of Farrel and clarke. Only the elements can do that: */ + Sgo_oceanaverage=femmodel->SealevelriseOceanAverage(Sgo); + + /*Sg is the sum of the eustatic term, and the ocean terms: */ + Sg_eustatic->Copy(Sg); Sg->AXPY(Sgo,1); + Sg->Shift(-Sgo_oceanaverage); + + /*convergence criterion:*/ + slrconvergence(&converged,Sg,Sg_old,eps_rel,eps_abs); + + /*Increase count: */ + count++; + if(converged==true){ + break; + } + if(count>=max_nonlinear_iterations){ + _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); + converged=true; + break; + } + + /*some minor verbosing adjustment:*/ + if(count>1)verboseconvolution=false; + + } + if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n"); + + xDelete(latitude); + xDelete(longitude); + xDelete(radius); + delete Sg_old; + + return Sg; +} + +void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ + + bool converged=true; + IssmDouble ndS,nS; + Vector *dSg = NULL; + + //compute norm(du) and norm(u) if requested + dSg=Sg_old->Duplicate(); Sg_old->Copy(dSg); dSg->AYPX(Sg,-1.0); + ndS=dSg->Norm(NORM_TWO); + + if(!xIsNan(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; + +} /*}}}*/