Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp (revision 22159) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp (revision 22160) @@ -58,6 +58,8 @@ Element* element=(Element*)elements->GetObjectByOffset(counter); element->InputUpdateFromConstant(0.0,SealevelriseCumDeltathicknessEnum); - element->InputUpdateFromConstant(0.0,SealevelRateEnum); - element->InputUpdateFromConstant(0.0,SealevelUmotionRateEnum); + element->InputUpdateFromConstant(0.0,SealevelNEsaRateEnum); + element->InputUpdateFromConstant(0.0,SealevelNGiaRateEnum); + element->InputUpdateFromConstant(0.0,SealevelUEsaRateEnum); + element->InputUpdateFromConstant(0.0,SealevelUGiaRateEnum); counter++; } Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.cpp (revision 22159) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.cpp (revision 22160) @@ -1641,13 +1641,20 @@ name==MaterialsRheologyEsbarEnum || name==SealevelEnum || - name==SealevelUmotionEnum || - name==SealevelUmotionRateEnum || - name==SealevelNmotionEnum || - name==SealevelEmotionEnum || - name==SealevelAbsoluteEnum || - name==SealevelEustaticEnum || + name==SealevelUEsaEnum || + name==SealevelUEsaRateEnum || + name==SealevelNEsaEnum || + name==SealevelNEsaRateEnum || + name==SealevelUNorthEsaEnum || + name==SealevelUEastEsaEnum || + name==SealevelRSLEustaticEnum || + name==SealevelRSLEustaticRateEnum || + name==SealevelRSLEnum || + name==SealevelRSLRateEnum || + name==SealevelUGiaEnum || + name==SealevelUGiaRateEnum || + name==SealevelNGiaEnum || + name==SealevelNGiaRateEnum || name==SealevelriseDeltathicknessEnum || name==SealevelriseCumDeltathicknessEnum || - name==SealevelRateEnum || name==EsaUmotionEnum || name==EsaNmotionEnum || Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/FemModel.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/FemModel.cpp (revision 22159) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/FemModel.cpp (revision 22160) @@ -2422,5 +2422,5 @@ #endif #ifdef _HAVE_SEALEVELRISE_ -void FemModel::SealevelriseEustatic(Vector* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/ +void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/ /*serialized vectors:*/ @@ -2461,13 +2461,13 @@ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea); + element->SealevelriseEustatic(pRSLgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea); eustatic_cpu+=eustatic_cpu_e; } - if(i%loop==0)pSgi->Assemble(); + if(i%loop==0)pRSLgi->Assemble(); } if(VerboseConvergence())_printf0_("\n"); /*One last time: */ - pSgi->Assemble(); + pRSLgi->Assemble(); /*Sum all eustatic components from all cpus:*/ @@ -2481,8 +2481,8 @@ } /*}}}*/ -void FemModel::SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/ +void FemModel::SealevelriseNonEustatic(Vector* pRSLgo, Vector* pRSLg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/ /*serialized vectors:*/ - IssmDouble* Sg_old=NULL; + IssmDouble* RSLg_old=NULL; IssmDouble eartharea=0; @@ -2492,5 +2492,5 @@ /*Serialize vectors from previous iteration:*/ - Sg_old=pSg_old->ToMPISerial(); + RSLg_old=pRSLg_old->ToMPISerial(); /*Go through elements, and add contribution from each element to the deflection vector wg:*/ @@ -2515,18 +2515,18 @@ if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf_("\r" << " convolution progress: " << (double)i/(double)ns*100 << "% "); Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseNonEustatic(pSgo,Sg_old,latitude,longitude,radius,eartharea); - } - if(i%loop==0)pSgo->Assemble(); + element->SealevelriseNonEustatic(pRSLgo,RSLg_old,latitude,longitude,radius,eartharea); + } + if(i%loop==0)pRSLgo->Assemble(); } if(verboseconvolution)if(VerboseConvergence())_printf_("\n"); /*Free ressources:*/ - xDelete(Sg_old); -} -/*}}}*/ -void FemModel::SealevelriseRotationalFeedback(Vector* pSgo_rot, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ + xDelete(RSLg_old); +} +/*}}}*/ +void FemModel::SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ /*serialized vectors:*/ - IssmDouble* Sg_old=NULL; + IssmDouble* RSLg_old=NULL; IssmDouble eartharea=0; IssmDouble eartharea_cpu=0; @@ -2537,5 +2537,5 @@ /*Serialize vectors from previous iteration:*/ - Sg_old=pSg_old->ToMPISerial(); + RSLg_old=pRSLg_old->ToMPISerial(); /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */ @@ -2551,5 +2551,5 @@ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseMomentOfInertia(&moi_list[0],Sg_old,eartharea); + element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,eartharea); moi_list_cpu[0] += moi_list[0]; moi_list_cpu[1] += moi_list[1]; @@ -2595,19 +2595,19 @@ (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi))); - pSgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times + pRSLgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times } /*Assemble mesh velocity*/ - pSgo_rot->Assemble(); + pRSLgo_rot->Assemble(); /*Free ressources:*/ - xDelete(Sg_old); - -} -/*}}}*/ -void FemModel::SealevelriseGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/ + xDelete(RSLg_old); + +} +/*}}}*/ +void FemModel::SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/ /*serialized vectors:*/ - IssmDouble* Sg=NULL; + IssmDouble* RSLg=NULL; IssmDouble eartharea=0; @@ -2617,5 +2617,5 @@ /*Serialize vectors from previous iteration:*/ - Sg=pSg->ToMPISerial(); + RSLg=pRSLg->ToMPISerial(); /*Go through elements, and add contribution from each element to the deflection vector wg:*/ @@ -2639,5 +2639,5 @@ if(VerboseConvergence())if(i%100==0)_printf0_("\r" << " convolution progress: " << (double)i/(double)ns*100 << "% "); Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz); + element->SealevelriseGeodetic(pUp,pNorth,pEast,RSLg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz); } if(i%loop==0){ @@ -2659,15 +2659,15 @@ /*Free ressources:*/ - xDelete(Sg); -} -/*}}}*/ -IssmDouble FemModel::SealevelriseOceanAverage(Vector* Sg) { /*{{{*/ - - IssmDouble* Sg_serial=NULL; + xDelete(RSLg); +} +/*}}}*/ +IssmDouble FemModel::SealevelriseOceanAverage(Vector* RSLg) { /*{{{*/ + + IssmDouble* RSLg_serial=NULL; IssmDouble oceanvalue,oceanvalue_cpu; IssmDouble oceanarea,oceanarea_cpu; /*Serialize vectors from previous iteration:*/ - Sg_serial=Sg->ToMPISerial(); + RSLg_serial=RSLg->ToMPISerial(); /*Initialize:*/ @@ -2679,5 +2679,5 @@ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); oceanarea_cpu += element->OceanArea(); - oceanvalue_cpu += element->OceanAverage(Sg_serial); + oceanvalue_cpu += element->OceanAverage(RSLg_serial); } ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); @@ -2688,5 +2688,5 @@ /*Free ressources:*/ - xDelete(Sg_serial); + xDelete(RSLg_serial); return oceanvalue/oceanarea; Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/FemModel.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/FemModel.h (revision 22159) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/FemModel.h (revision 22160) @@ -129,5 +129,5 @@ void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop); void SealevelriseRotationalFeedback(Vector* pSgo_rot, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); - void SealevelriseGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); + void SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); IssmDouble SealevelriseOceanAverage(Vector* Sg); #endif Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h (revision 22159) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h (revision 22160) @@ -54,6 +54,8 @@ void steric_core(FemModel* femmodel); Vector* sealevelrise_core_eustatic(FemModel* femmodel); -Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* Sg_eustatic); -void sealevelrise_core_elastic(Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* Sg); +Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* RSLg_eustatic); +void sealevelrise_core_elastic(Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* RSLg); +void sealevelrise_core_viscous(Vector** pU_gia,Vector** pN_gia,FemModel* femmodel,Vector* RSLg); +void sealevelrise_diagnostics(FemModel* femmodel,Vector* RSLg); IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel); @@ -69,5 +71,5 @@ void TransferSealevel(FemModel* femmodel,int forcingenum); void EarthMassTransport(FemModel* femmodel); -void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs); +void slrconvergence(bool* pconverged, Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs); //solution configuration Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp (revision 22159) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp (revision 22160) @@ -61,15 +61,22 @@ /*variables:*/ - Vector *Sg = NULL; - Vector *Sg_rate = NULL; - Vector *Sg_eustatic = NULL; - Vector *U_radial = NULL; - Vector *U_radial_rate = NULL; - Vector *U_north = NULL; - Vector *U_east = NULL; - Vector* cumdeltathickness=NULL; + Vector *RSLg = NULL; + Vector *RSLg_rate = NULL; + Vector *RSLg_eustatic = NULL; + Vector *RSLg_eustatic_rate = NULL; + Vector *U_esa = NULL; + Vector *U_esa_rate = NULL; + Vector *N_esa = NULL; + Vector *N_esa_rate = NULL; + Vector *U_north_esa = NULL; + Vector *U_east_esa = NULL; + Vector *N_gia= NULL; + Vector *U_gia= NULL; + Vector *N_gia_rate= NULL; + Vector *U_gia_rate= NULL; /*parameters:*/ bool iscoupler; + int solution_type; int configuration_type; int modelid,earthid; @@ -93,4 +100,5 @@ femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum); femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); if(iscoupler){ @@ -132,25 +140,50 @@ /*call eustatic core (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */ - Sg_eustatic=sealevelrise_core_eustatic(femmodel); + RSLg_eustatic=sealevelrise_core_eustatic(femmodel); /*call non-eustatic core (ocean loading tems - 2nd and 5th terms on the RHS of Farrel and Clark) */ - Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); - + RSLg=sealevelrise_core_noneustatic(femmodel,RSLg_eustatic); /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */ - sealevelrise_core_elastic(&U_radial,&U_north,&U_east,femmodel,Sg); + sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg); + + /*compute viscosus (GIA) geodetic signatures:*/ + sealevelrise_core_viscous(&U_gia,&N_gia,femmodel,RSLg); + + /*compute sea-level rise (low-order spherical harmonics coefficients) diagnostics:*/ + sealevelrise_diagnostics(femmodel,RSLg); + + /*recover N_esa = U_esa + RSLg:*/ + N_esa=U_esa->Duplicate(); U_esa->Copy(N_esa); N_esa->AXPY(RSLg,1); /*transform these values into rates (as we only run this once each frequency turn:*/ - Sg_rate=Sg->Duplicate(); Sg->Copy(Sg_rate); Sg_rate->Scale(1/(dt*frequency)); - U_radial_rate=U_radial->Duplicate(); U_radial->Copy(U_radial_rate); U_radial_rate->Scale(1/(dt*frequency)); - - - /*get some results into elements:*/ - InputUpdateFromVectorx(femmodel,Sg_rate,SealevelRateEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,U_radial_rate,SealevelUmotionRateEnum,VertexSIdEnum); - if (horiz){ - InputUpdateFromVectorx(femmodel,U_north,SealevelNmotionEnum,VertexSIdEnum); // north motion - InputUpdateFromVectorx(femmodel,U_east,SealevelEmotionEnum,VertexSIdEnum); // east motion - } + N_esa_rate=N_esa->Duplicate(); N_esa->Copy(N_esa_rate); N_esa_rate->Scale(1/(dt*frequency)); + U_esa_rate=U_esa->Duplicate(); U_esa->Copy(U_esa_rate); U_esa_rate->Scale(1/(dt*frequency)); + N_gia_rate=N_gia->Duplicate(); N_gia->Copy(N_gia_rate); N_gia_rate->Scale(1/(dt*frequency)); + U_gia_rate=U_gia->Duplicate(); U_gia->Copy(U_gia_rate); U_gia_rate->Scale(1/(dt*frequency)); + RSLg_rate=RSLg->Duplicate(); RSLg->Copy(RSLg_rate); RSLg_rate->Scale(1/(dt*frequency)); + RSLg_eustatic_rate=RSLg_eustatic->Duplicate(); RSLg_eustatic->Copy(RSLg_eustatic_rate); RSLg_eustatic_rate->Scale(1/(dt*frequency)); + + /*get some results into elements:{{{*/ + if(solution_type==SealevelriseSolutionEnum){ + InputUpdateFromVectorx(femmodel,U_esa,SealevelUEsaEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,N_esa,SealevelNEsaEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,U_gia,SealevelUGiaEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,N_gia,SealevelNGiaEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,RSLg,SealevelRSLEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,RSLg_eustatic,SealevelRSLEustaticEnum,VertexSIdEnum); + if (horiz){ + InputUpdateFromVectorx(femmodel,U_north_esa,SealevelUNorthEsaEnum,VertexSIdEnum); // north motion + InputUpdateFromVectorx(femmodel,U_east_esa,SealevelUEsaEnum,VertexSIdEnum); // east motion + } + } + else{ + InputUpdateFromVectorx(femmodel,U_esa_rate,SealevelUEsaRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,N_esa_rate,SealevelNEsaRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,U_gia_rate,SealevelUGiaRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,N_gia_rate,SealevelNGiaRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,RSLg_rate,SealevelRSLRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,RSLg_eustatic_rate,SealevelRSLEustaticRateEnum,VertexSIdEnum); + } /*}}}*/ } @@ -158,6 +191,8 @@ if(iscoupler){ /*transfer sea level back to ice caps:*/ - TransferSealevel(femmodel,SealevelRateEnum); - TransferSealevel(femmodel,SealevelUmotionRateEnum); + TransferSealevel(femmodel,SealevelNEsaRateEnum); + TransferSealevel(femmodel,SealevelNGiaRateEnum); + TransferSealevel(femmodel,SealevelUEsaRateEnum); + TransferSealevel(femmodel,SealevelUGiaRateEnum); //reset cumdeltathickness to 0: @@ -168,14 +203,24 @@ femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter. - /*Free ressources:*/ - delete Sg_eustatic; - delete Sg; - delete Sg_rate; - delete U_radial; - delete U_radial_rate; + /*free ressources:{{{*/ + delete RSLg; + delete RSLg_rate; + delete RSLg_eustatic; + delete RSLg_eustatic_rate; + delete U_esa; + delete U_esa_rate; + delete N_esa; + delete N_esa_rate; + if(horiz){ - delete U_north; - delete U_east; - } + delete U_north_esa; + delete U_east_esa; + } + delete N_gia; + delete U_gia; + delete N_gia_rate; + delete U_gia_rate; + /*}}}*/ + } /*}}}*/ @@ -184,8 +229,10 @@ /*variables:*/ Vector *bedrock = NULL; - Vector *Sg_absolute = NULL; + Vector *SL = NULL; Vector *steric_rate_g = NULL; - Vector *rsl_g = NULL; - Vector *bedrock_rate_g = NULL; + Vector *U_esa_rate= NULL; + Vector *N_esa_rate= NULL; + Vector *U_gia_rate= NULL; + Vector *N_gia_rate= NULL; /*parameters: */ @@ -210,33 +257,44 @@ if(VerboseSolution()) _printf0_(" computing steric sea level rise\n"); - /*Retrieve absolute sea level, RSL rate, bedrock uplift rate and steric rate, as vectors:*/ + /*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/ GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum); - GetVectorFromInputsx(&Sg_absolute,femmodel,SealevelEnum,VertexSIdEnum); + GetVectorFromInputsx(&SL,femmodel,SealevelEnum,VertexSIdEnum); GetVectorFromInputsx(&steric_rate_g,femmodel,SealevelriseStericRateEnum,VertexSIdEnum); if(geodetic){ - GetVectorFromInputsx(&rsl_g,femmodel,SealevelRateEnum,VertexSIdEnum); - GetVectorFromInputsx(&bedrock_rate_g,femmodel,SealevelUmotionRateEnum,VertexSIdEnum); - } - - /*compute: absolute sea level change = initial absolute sea level + relative sea level change rate * dt + vertical motion rate * dt + steric_rate * dt*/ - if(rsl_g) Sg_absolute->AXPY(rsl_g,dt); - if(bedrock_rate_g)Sg_absolute->AXPY(bedrock_rate_g,dt); - Sg_absolute->AXPY(steric_rate_g,dt); + GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&U_gia_rate,femmodel,SealevelUGiaRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&N_esa_rate,femmodel,SealevelNEsaRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&N_gia_rate,femmodel,SealevelNGiaRateEnum,VertexSIdEnum); + } + + /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate - U_gia_rate - U_esa_rate ) * dt + steric_rate * dt*/ + if(geodetic){ + SL->AXPY(N_gia_rate,dt); + SL->AXPY(N_esa_rate,dt); + SL->AXPY(U_gia_rate,-dt); //notice -dt! + SL->AXPY(U_esa_rate,-dt); //same here + } + SL->AXPY(steric_rate_g,dt); /*compute new bedrock position: */ - if(bedrock_rate_g)bedrock->AXPY(bedrock_rate_g,dt); + if(geodetic){ + bedrock->AXPY(U_esa_rate,dt); + bedrock->AXPY(U_gia_rate,dt); + } /*update element inputs:*/ InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,bedrock,SealevelUmotionEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,SL,SealevelEnum,VertexSIdEnum); /*Free ressources:*/ delete bedrock; - delete Sg_absolute; + delete SL; delete steric_rate_g; - if(rsl_g)delete rsl_g; - if(bedrock_rate_g)delete bedrock_rate_g; - + if(geodetic){ + delete U_esa_rate; + delete U_gia_rate; + delete N_esa_rate; + delete N_gia_rate; + } } /*}}}*/ @@ -245,6 +303,6 @@ /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/ - Vector *Sgi = NULL; - IssmDouble Sgi_oceanaverage = 0; + Vector *RSLgi = NULL; + IssmDouble RSLgi_oceanaverage = 0; /*parameters: */ @@ -271,18 +329,15 @@ /*Initialize:*/ - Sgi = new Vector(gsize); + RSLgi = new Vector(gsize); /*call the eustatic main module: */ - femmodel->SealevelriseEustatic(Sgi,&eustatic, latitude, longitude, radius,loop); //this computes - - /*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 + femmodel->SealevelriseEustatic(RSLgi,&eustatic, latitude, longitude, radius,loop); //this computes + + /*we need to average RSLgi over the ocean: RHS term 4 in Eq.4 of Farrel and clarke. Only the elements can do that: */ + RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi); + + /*RSLg 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); - - /*save eustatic value for results: */ - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelEustaticEnum,-eustatic)); + RSLgi->Shift(-eustatic-RSLgi_oceanaverage); /*clean up and return:*/ @@ -290,7 +345,7 @@ xDelete(longitude); xDelete(radius); - return Sgi; + return RSLgi; }/*}}}*/ -Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* Sg_eustatic){ /*{{{*/ +Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* RSLg_eustatic){ /*{{{*/ /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5. @@ -298,10 +353,10 @@ - Vector *Sg = NULL; - Vector *Sg_old = NULL; - - Vector *Sgo = NULL; //ocean convolution of the perturbation to gravity potential. - Vector *Sgo_rot= NULL; // rotational feedback - IssmDouble Sgo_oceanaverage = 0; //average of Sgo over the ocean. + Vector *RSLg = NULL; + Vector *RSLg_old = NULL; + + Vector *RSLgo = NULL; //ocean convolution of the perturbation to gravity potential. + Vector *RSLgo_rot= NULL; // rotational feedback + IssmDouble RSLgo_oceanaverage = 0; //average of RSLgo over the ocean. /*parameters: */ @@ -340,10 +395,10 @@ /*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(); + RSLg = new Vector(gsize); + RSLg->Assemble(); + RSLg_eustatic->Copy(RSLg); //first initialize RSLg with the eustatic component computed in sealevelrise_core_eustatic. + + RSLg_old = new Vector(gsize); + RSLg_old->Assemble(); count=1; @@ -354,37 +409,37 @@ //save pointer to old sea level rise - delete Sg_old; Sg_old=Sg; + delete RSLg_old; RSLg_old=RSLg; /*Initialize solution vector: */ - Sg = new Vector(gsize); Sg->Assemble(); - Sgo = new Vector(gsize); Sgo->Assemble(); + RSLg = new Vector(gsize); RSLg->Assemble(); + RSLgo = new Vector(gsize); RSLgo->Assemble(); /*call the non eustatic module: */ - femmodel->SealevelriseNonEustatic(Sgo, Sg_old, latitude, longitude, radius,verboseconvolution,loop); + femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, latitude, longitude, radius,verboseconvolution,loop); /*assemble solution vector: */ - Sgo->Assemble(); + RSLgo->Assemble(); if(rotation){ /*call rotational feedback module: */ - Sgo_rot = new Vector(gsize); Sgo_rot->Assemble(); - femmodel->SealevelriseRotationalFeedback(Sgo_rot,Sg_old,latitude,longitude,radius); - Sgo_rot->Assemble(); - - Sgo->AXPY(Sgo_rot,1); - } - - /*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); + RSLgo_rot = new Vector(gsize); RSLgo_rot->Assemble(); + femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,latitude,longitude,radius); + RSLgo_rot->Assemble(); + + RSLgo->AXPY(RSLgo_rot,1); + } + + /*we need to average RSLgo over the ocean: RHS term 5 in Eq.4 of Farrel and clarke. Only the elements can do that: */ + RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo); + + /*RSLg is the sum of the eustatic term, and the ocean terms: */ + RSLg_eustatic->Copy(RSLg); RSLg->AXPY(RSLgo,1); + RSLg->Shift(-RSLgo_oceanaverage); /*convergence criterion:*/ - slrconvergence(&converged,Sg,Sg_old,eps_rel,eps_abs); + slrconvergence(&converged,RSLg,RSLg_old,eps_rel,eps_abs); /*free ressources: */ - delete Sgo; + delete RSLgo; /*Increase count: */ @@ -408,13 +463,13 @@ xDelete(longitude); xDelete(radius); - delete Sg_old; - - return Sg; + delete RSLg_old; + + return RSLg; } /*}}}*/ -void sealevelrise_core_elastic(Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* Sg){ /*{{{*/ - - Vector *U_radial = NULL; - Vector *U_north = NULL; - Vector *U_east = NULL; +void sealevelrise_core_elastic(Vector** pU_esa, Vector** pU_north_esa,Vector** pU_east_esa,FemModel* femmodel,Vector* RSLg){ /*{{{*/ + + Vector *U_esa = NULL; + Vector *U_north_esa = NULL; + Vector *U_east_esa = NULL; /*parameters: */ @@ -441,8 +496,8 @@ /*intialize vectors:*/ - U_radial = new Vector(gsize); + U_esa = new Vector(gsize); if (horiz){ - U_north = new Vector(gsize); - U_east = new Vector(gsize); + U_north_esa = new Vector(gsize); + U_east_esa = new Vector(gsize); } @@ -452,11 +507,11 @@ /*call the elastic main modlule:*/ - femmodel->SealevelriseGeodetic(U_radial,U_north,U_east,Sg,latitude,longitude,radius,xx,yy,zz,loop,horiz); + femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,latitude,longitude,radius,xx,yy,zz,loop,horiz); /*Assign output pointers:*/ - *pU_radial=U_radial; + *pU_esa=U_esa; if(horiz){ - *pU_east=U_east; - *pU_north=U_north; + *pU_east_esa=U_east_esa; + *pU_north_esa=U_north_esa; } @@ -468,4 +523,34 @@ delete yy; delete zz; +} +/*}}}*/ +void sealevelrise_core_viscous(Vector** pU_gia, Vector** pN_gia,FemModel* femmodel,Vector* RSLg){ /*{{{*/ + + Vector *U_gia = NULL; + Vector *N_gia = NULL; + + int configuration_type; + int gsize; + + /*retrieve some parameters:*/ + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + + /*find size of vectors:*/ + gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); + + /*intialize vectors:*/ + U_gia = new Vector(gsize); U_gia->Set(0); + N_gia = new Vector(gsize); N_gia->Set(0); + + /*Assign output pointers:*/ + *pU_gia=U_gia; + *pN_gia=N_gia; + +} +/*}}}*/ +void sealevelrise_diagnostics(FemModel* femmodel,Vector* RSLg){ /*{{{*/ + + /*compute spherical harmonics deg 1 and deg 2 coefficeints:*/ + } /*}}}*/ @@ -748,18 +833,18 @@ } /*}}}*/ -void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ +void slrconvergence(bool* pconverged, Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ bool converged=true; IssmDouble ndS,nS; - Vector *dSg = NULL; + Vector *dRSLg = 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); + dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0); + ndS=dRSLg->Norm(NORM_TWO); if (xIsNan(ndS)) _error_("convergence criterion is NaN!"); if(!xIsNan(eps_rel)){ - nS=Sg_old->Norm(NORM_TWO); + nS=RSLg_old->Norm(NORM_TWO); if (xIsNan(nS)) _error_("convergence criterion is NaN!"); } @@ -767,5 +852,5 @@ //clean up - delete dSg; + delete dRSLg; //print Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumDefinitions.h (revision 22159) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumDefinitions.h (revision 22160) @@ -783,11 +783,18 @@ /*Sea Level Rise{{{*/ SealevelEnum, - SealevelRateEnum, - SealevelUmotionEnum, - SealevelUmotionRateEnum, - SealevelNmotionEnum, - SealevelEmotionEnum, - SealevelAbsoluteEnum, - SealevelEustaticEnum, + SealevelUEsaEnum, + SealevelUEsaRateEnum, + SealevelNEsaEnum, + SealevelNEsaRateEnum, + SealevelUNorthEsaEnum, + SealevelUEastEsaEnum, + SealevelUGiaEnum, + SealevelUGiaRateEnum, + SealevelNGiaEnum, + SealevelNGiaRateEnum, + SealevelRSLEustaticEnum, + SealevelRSLEustaticRateEnum, + SealevelRSLEnum, + SealevelRSLRateEnum, SealevelObsEnum, SealevelWeightsEnum, Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumToStringx.cpp (revision 22159) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/EnumToStringx.cpp (revision 22160) @@ -763,11 +763,18 @@ case LevelsetReinitFrequencyEnum : return "LevelsetReinitFrequency"; case SealevelEnum : return "Sealevel"; - case SealevelRateEnum : return "SealevelRate"; - case SealevelUmotionEnum : return "SealevelUmotion"; - case SealevelUmotionRateEnum : return "SealevelUmotionRate"; - case SealevelNmotionEnum : return "SealevelNmotion"; - case SealevelEmotionEnum : return "SealevelEmotion"; - case SealevelAbsoluteEnum : return "SealevelAbsolute"; - case SealevelEustaticEnum : return "SealevelEustatic"; + case SealevelUEsaEnum : return "SealevelUEsa"; + case SealevelUEsaRateEnum : return "SealevelUEsaRate"; + case SealevelNEsaEnum : return "SealevelNEsa"; + case SealevelNEsaRateEnum : return "SealevelNEsaRate"; + case SealevelUNorthEsaEnum : return "SealevelUNorthEsa"; + case SealevelUEastEsaEnum : return "SealevelUEastEsa"; + case SealevelUGiaEnum : return "SealevelUGia"; + case SealevelUGiaRateEnum : return "SealevelUGiaRate"; + case SealevelNGiaEnum : return "SealevelNGia"; + case SealevelNGiaRateEnum : return "SealevelNGiaRate"; + case SealevelRSLEustaticEnum : return "SealevelRSLEustatic"; + case SealevelRSLEustaticRateEnum : return "SealevelRSLEustaticRate"; + case SealevelRSLEnum : return "SealevelRSL"; + case SealevelRSLRateEnum : return "SealevelRSLRate"; case SealevelObsEnum : return "SealevelObs"; case SealevelWeightsEnum : return "SealevelWeights"; Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/StringToEnumx.cpp (revision 22159) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/shared/Enum/StringToEnumx.cpp (revision 22160) @@ -781,11 +781,18 @@ else if (strcmp(name,"LevelsetReinitFrequency")==0) return LevelsetReinitFrequencyEnum; else if (strcmp(name,"Sealevel")==0) return SealevelEnum; - else if (strcmp(name,"SealevelRate")==0) return SealevelRateEnum; - else if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum; - else if (strcmp(name,"SealevelUmotionRate")==0) return SealevelUmotionRateEnum; - 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,"SealevelUEsa")==0) return SealevelUEsaEnum; + else if (strcmp(name,"SealevelUEsaRate")==0) return SealevelUEsaRateEnum; + else if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum; + else if (strcmp(name,"SealevelNEsaRate")==0) return SealevelNEsaRateEnum; + else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum; + else if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum; + else if (strcmp(name,"SealevelUGia")==0) return SealevelUGiaEnum; + else if (strcmp(name,"SealevelUGiaRate")==0) return SealevelUGiaRateEnum; + else if (strcmp(name,"SealevelNGia")==0) return SealevelNGiaEnum; + else if (strcmp(name,"SealevelNGiaRate")==0) return SealevelNGiaRateEnum; + else if (strcmp(name,"SealevelRSLEustatic")==0) return SealevelRSLEustaticEnum; + else if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum; + else if (strcmp(name,"SealevelRSL")==0) return SealevelRSLEnum; + else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum; else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum; else if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum; @@ -868,5 +875,8 @@ else if (strcmp(name,"Node")==0) return NodeEnum; else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; - else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum; else if (strcmp(name,"Param")==0) return ParamEnum; else if (strcmp(name,"Moulin")==0) return MoulinEnum; @@ -875,8 +885,5 @@ else if (strcmp(name,"Profiler")==0) return ProfilerEnum; else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"Masscon")==0) return MassconEnum; + else if (strcmp(name,"Masscon")==0) return MassconEnum; else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; @@ -991,5 +998,8 @@ else if (strcmp(name,"IceMass")==0) return IceMassEnum; else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; - else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum; else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum; else if (strcmp(name,"TotalGroundedBmb")==0) return TotalGroundedBmbEnum; @@ -998,8 +1008,5 @@ else if (strcmp(name,"P0Array")==0) return P0ArrayEnum; else if (strcmp(name,"P1")==0) return P1Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"P1DG")==0) return P1DGEnum; + else if (strcmp(name,"P1DG")==0) return P1DGEnum; else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum; else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;