Index: /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 24983) +++ /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 24984) @@ -24,9 +24,4 @@ if(world_rank!=0){ femmodel_init= new FemModel(argc,argv,evaluation_comm); - - /*hack: we don't have any place to put solution modules that fit between the model processor - *and the cores:*/ - //for (int i=0;inummodels;i++) if(femmodel_init->analysis_type_list[i]==SealevelriseAnalysisEnum)sealevelrise_core_geometry(femmodel_init); - femmodel_init->profiler->Start(CORE); } Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24983) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24984) @@ -4673,14 +4673,4 @@ #endif #ifdef _HAVE_SEALEVELRISE_ -void FemModel::SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz) { /*{{{*/ - - /*Run sealevelrie geometry routine in elements:*/ - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseGeometry(latitude,longitude,radius,xx,yy,zz); - } - -} -/*}}}*/ void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks) { /*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24983) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24984) @@ -164,5 +164,4 @@ #ifdef _HAVE_SEALEVELRISE_ void SealevelriseEustatic(Vector* pSgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks); - void SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, SealevelMasks* masks,bool verboseconvolution); void SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks); Index: /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 24983) +++ /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 24984) @@ -60,5 +60,5 @@ break; case SealevelriseSolutionEnum: - solutioncore=&sealevelrise_core; + solutioncore=&sealevelchange_core; break; case EsaSolutionEnum: Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 24983) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 24984) @@ -53,7 +53,7 @@ void bmb_core(FemModel* femmodel); void damage_core(FemModel* femmodel); -void sealevelrise_core(FemModel* femmodel); -void geodetic_core(FemModel* femmodel); -void steric_core(FemModel* femmodel); +void sealevelchange_core(FemModel* femmodel); +void grd_core(FemModel* femmodel); +void dynstr_core(FemModel* femmodel); void sealevelrise_core_geometry(FemModel* femmodel); SealevelMasks* sealevelrise_core_masks(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24983) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24984) @@ -1,4 +1,4 @@ -/*!\file: sealevelrise_core.cpp - * \brief: core of the SLR solution +/*!\file: sealevelchange_core.cpp + * \brief: core of the sea-level change solution */ @@ -14,6 +14,6 @@ -/*cores:*/ -void sealevelrise_core(FemModel* femmodel){ /*{{{*/ +/*main cores:*/ +void sealevelchange_core(FemModel* femmodel){ /*{{{*/ /*Start profiler*/ @@ -46,8 +46,8 @@ /*Run geodetic:*/ - geodetic_core(femmodel); + grd_core(femmodel); /*Run steric core for sure:*/ - steric_core(femmodel); + dynstr_core(femmodel); /*Save results: */ @@ -69,5 +69,7 @@ } /*}}}*/ -void geodetic_core(FemModel* femmodel){ /*{{{*/ +void grd_core(FemModel* femmodel){ /*{{{*/ + + /*Gravity rotation deformation core GRD: */ /*variables:*/ @@ -242,4 +244,85 @@ } /*}}}*/ +void dynstr_core(FemModel* femmodel){ /*{{{*/ + + /*variables:*/ + Vector *bedrock = NULL; + Vector *SL = NULL; + Vector *steric_rate_g = NULL; + Vector *dynamic_rate_g = NULL; + Vector *hydro_rate_g = NULL; + Vector *U_esa_rate= NULL; + Vector *N_esa_rate= NULL; + Vector *U_gia_rate= NULL; + Vector *N_gia_rate= NULL; + + /*parameters: */ + bool isslr=0; + int solution_type; + IssmDouble dt; + int geodetic=0; + + /*Retrieve parameters:*/ + femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); + + /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/ + if(solution_type==SealevelriseSolutionEnum)isslr=1; + + /*Should we be here?:*/ + if(!isslr)return; + + /*Verbose: */ + if(VerboseSolution()) _printf0_(" computing steric sea level rise\n"); + + /*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/ + GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum); + GetVectorFromInputsx(&SL,femmodel,SealevelEnum,VertexSIdEnum); + GetStericRate(&steric_rate_g,femmodel); + GetDynamicRate(&dynamic_rate_g,femmodel); + GetVectorFromInputsx(&hydro_rate_g,femmodel,SealevelriseHydroRateEnum,VertexSIdEnum); + if(geodetic){ + 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) * dt + steric_rate + dynamic_rate + hydro_rate* dt*/ + if(geodetic){ + SL->AXPY(N_gia_rate,dt); + SL->AXPY(N_esa_rate,dt); + } + SL->AXPY(steric_rate_g,dt); + SL->AXPY(dynamic_rate_g,dt); + SL->AXPY(hydro_rate_g,dt); + + /*compute new bedrock position: */ + if(geodetic){ + bedrock->AXPY(U_esa_rate,dt); + bedrock->AXPY(U_gia_rate,dt); + } + + /*update element inputs:*/ + InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,SL,SealevelEnum,VertexSIdEnum); + + /*Free ressources:*/ + delete bedrock; + delete SL; + delete steric_rate_g; + delete dynamic_rate_g; + delete hydro_rate_g; + if(geodetic){ + delete U_esa_rate; + delete U_gia_rate; + delete N_esa_rate; + delete N_gia_rate; + } +} +/*}}}*/ + SealevelMasks* sealevelrise_core_masks(FemModel* femmodel) { /*{{{*/ @@ -289,6 +372,10 @@ if(horiz) VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); - /*call the FemModel geometry module: */ - femmodel->SealevelriseGeometry(latitude, longitude, xx,yy,zz,radius); + + /*Run sealevelrie geometry routine in elements:*/ + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->SealevelriseGeometry(latitude,longitude,radius,xx,yy,zz); + } /*Free ressources:*/ @@ -652,84 +739,4 @@ } /*}}}*/ -void steric_core(FemModel* femmodel){ /*{{{*/ - - /*variables:*/ - Vector *bedrock = NULL; - Vector *SL = NULL; - Vector *steric_rate_g = NULL; - Vector *dynamic_rate_g = NULL; - Vector *hydro_rate_g = NULL; - Vector *U_esa_rate= NULL; - Vector *N_esa_rate= NULL; - Vector *U_gia_rate= NULL; - Vector *N_gia_rate= NULL; - - /*parameters: */ - bool isslr=0; - int solution_type; - IssmDouble dt; - int geodetic=0; - - /*Retrieve parameters:*/ - femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); - femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); - femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); - - /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/ - if(solution_type==SealevelriseSolutionEnum)isslr=1; - - /*Should we be here?:*/ - if(!isslr)return; - - /*Verbose: */ - if(VerboseSolution()) _printf0_(" computing steric sea level rise\n"); - - /*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/ - GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum); - GetVectorFromInputsx(&SL,femmodel,SealevelEnum,VertexSIdEnum); - GetStericRate(&steric_rate_g,femmodel); - GetDynamicRate(&dynamic_rate_g,femmodel); - GetVectorFromInputsx(&hydro_rate_g,femmodel,SealevelriseHydroRateEnum,VertexSIdEnum); - if(geodetic){ - 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) * dt + steric_rate + dynamic_rate + hydro_rate* dt*/ - if(geodetic){ - SL->AXPY(N_gia_rate,dt); - SL->AXPY(N_esa_rate,dt); - } - SL->AXPY(steric_rate_g,dt); - SL->AXPY(dynamic_rate_g,dt); - SL->AXPY(hydro_rate_g,dt); - - /*compute new bedrock position: */ - if(geodetic){ - bedrock->AXPY(U_esa_rate,dt); - bedrock->AXPY(U_gia_rate,dt); - } - - /*update element inputs:*/ - InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum); - InputUpdateFromVectorx(femmodel,SL,SealevelEnum,VertexSIdEnum); - - /*Free ressources:*/ - delete bedrock; - delete SL; - delete steric_rate_g; - delete dynamic_rate_g; - delete hydro_rate_g; - if(geodetic){ - delete U_esa_rate; - delete U_gia_rate; - delete N_esa_rate; - delete N_gia_rate; - } -} -/*}}}*/ /*support routines:*/ Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 24983) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 24984) @@ -228,5 +228,5 @@ /*Sea level rise: */ - if(isslr) sealevelrise_core(femmodel); + if(isslr) sealevelchange_core(femmodel); /*unload results*/