Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 24920) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 24921) @@ -282,5 +282,6 @@ int dslmodel=0; - bool elastic=false; + IssmDouble* G_rigid = NULL; + IssmDouble* G_rigid_local = NULL; IssmDouble* G_elastic = NULL; IssmDouble* G_elastic_local = NULL; @@ -338,125 +339,132 @@ } /*}}}*/ /*Deal with elasticity {{{*/ - iomodel->FetchData(&elastic,"md.slr.elastic"); - if(elastic){ - - /*love numbers: */ - 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*/ - iomodel->FetchData(°acc,"md.slr.degacc"); - M=reCast(180./degacc+1.); - - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. - #ifdef _HAVE_AD_ - G_elastic=xNew(M,"t"); - U_elastic=xNew(M,"t"); - H_elastic=xNew(M,"t"); - #else - G_elastic=xNew(M); - U_elastic=xNew(M); - H_elastic=xNew(M); - #endif - - /*compute combined legendre + love number (elastic green function:*/ - m=DetermineLocalSize(M,IssmComm::GetComm()); - GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. - #ifdef _HAVE_AD_ - G_elastic_local=xNew(m,"t"); - U_elastic_local=xNew(m,"t"); - H_elastic_local=xNew(m,"t"); - #else - G_elastic_local=xNew(m); - U_elastic_local=xNew(m); - H_elastic_local=xNew(m); - #endif - - for(int i=lower_row;i(i)*degacc * PI / 180.0; - - G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0); - U_elastic_local[i-lower_row]= (love_h[nl-1])/2.0/sin(alpha/2.0); - H_elastic_local[i-lower_row]= 0; - IssmDouble Pn = 0.; - IssmDouble Pn1 = 0.; - IssmDouble Pn2 = 0.; - IssmDouble Pn_p = 0.; - IssmDouble Pn_p1 = 0.; - IssmDouble Pn_p2 = 0.; - - for (int n=0;n(IssmComm::GetSize()); - int* displs=xNew(IssmComm::GetSize()); - - //recvcounts: - ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); - - /*displs: */ - ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); - - /*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); - xDelete(displs); - - /*}}}*/ - - /*Avoid singularity at 0: */ - 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); - } /*}}}*/ + /*love numbers: */ + 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*/ + iomodel->FetchData(°acc,"md.slr.degacc"); + M=reCast(180./degacc+1.); + + // AD performance is sensitive to calls to ensurecontiguous. + // // Providing "t" will cause ensurecontiguous to be called. + #ifdef _HAVE_AD_ + G_rigid=xNew(M,"t"); + G_elastic=xNew(M,"t"); + U_elastic=xNew(M,"t"); + H_elastic=xNew(M,"t"); + #else + G_rigid=xNew(M); + G_elastic=xNew(M); + U_elastic=xNew(M); + H_elastic=xNew(M); + #endif + + /*compute combined legendre + love number (elastic green function:*/ + m=DetermineLocalSize(M,IssmComm::GetComm()); + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); + // AD performance is sensitive to calls to ensurecontiguous. + // // Providing "t" will cause ensurecontiguous to be called. + #ifdef _HAVE_AD_ + G_elastic_local=xNew(m,"t"); + G_rigid_local=xNew(m,"t"); + U_elastic_local=xNew(m,"t"); + H_elastic_local=xNew(m,"t"); + #else + G_elastic_local=xNew(m); + G_rigid_local=xNew(m); + U_elastic_local=xNew(m); + H_elastic_local=xNew(m); + #endif + + for(int i=lower_row;i(i)*degacc * PI / 180.0; + + G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); + G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row]; + U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row]; + H_elastic_local[i-lower_row]= 0; + IssmDouble Pn = 0.; + IssmDouble Pn1 = 0.; + IssmDouble Pn2 = 0.; + IssmDouble Pn_p = 0.; + IssmDouble Pn_p1 = 0.; + IssmDouble Pn_p2 = 0.; + + for (int n=0;n(IssmComm::GetSize()); + int* displs=xNew(IssmComm::GetSize()); + + //recvcounts: + ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); + + /*displs: */ + ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); + + /*All gather:*/ + ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + 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); + xDelete(displs); + + /*}}}*/ + + /*Avoid singularity at 0: */ + G_rigid[0]=G_rigid[1]; + parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); + 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_rigid); + xDelete(G_rigid_local); + xDelete(G_elastic); + xDelete(G_elastic_local); + xDelete(U_elastic); + xDelete(U_elastic_local); + xDelete(H_elastic); + xDelete(H_elastic_local); +/*}}}*/ /*Transitions:{{{ */ iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.slr.transitions"); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24920) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24921) @@ -5598,5 +5598,7 @@ /*elastic green function:*/ + int index; IssmDouble* G_elastic_precomputed=NULL; + IssmDouble* G_rigid_precomputed=NULL; int M; @@ -5649,7 +5651,9 @@ /*recover elastic green function:*/ + DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M); + if(computeelastic){ - DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); - _assert_(parameter); + parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M); } @@ -5748,25 +5752,17 @@ for(int i=0;i(alpha/PI*reCast(M-1)); //Elastic component (from Eq 17 in Adhikari et al, GMD 2015) - if(computeelastic){ - int index=reCast(alpha/PI*reCast(M-1)); - G_elastic += G_elastic_precomputed[index]; - } + if(computeelastic) G_elastic += G_elastic_precomputed[index]; /*Add all components to the Sgi or Sgo solution vectors:*/ - Sgi[i]+=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic); - + Sgi[i]+=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid_precomputed[index]+G_elastic); } } @@ -5940,10 +5936,10 @@ /*precomputed elastic green functions:*/ + IssmDouble* G_rigid_precomputed = NULL; IssmDouble* G_elastic_precomputed = NULL; int M; - - /*computation of Green functions:*/ - IssmDouble* G_elastic= NULL; - IssmDouble* G_rigid= NULL; + int index; + IssmDouble alpha; + IssmDouble delPhi,delLambda; /*optimization:*/ @@ -6020,46 +6016,34 @@ /*}}}*/ + /*recover rigid and elastic green functions:*/ + DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M); + if(computeelastic){ - /*recover elastic green function:*/ - DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); + parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M); - - /*initialize G_elastic:*/ - G_elastic=xNewZeroInit(gsize); - } - if(computerigid) G_rigid=xNewZeroInit(gsize); - - IssmDouble alpha; - IssmDouble delPhi,delLambda; - + } + + for(int i=0;i(alpha/PI*(M-1)); /*Rigid earth gravitational perturbation: */ if(computerigid){ - G_rigid[i]=1.0/2.0/sin(alpha/2.0); - //values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_rigid[i]; - Sgo[i]+=3*rho_water/rho_earth*area/eartharea*S*G_rigid[i]; + Sgo[i]+=3*rho_water/rho_earth*area/eartharea*S*G_rigid_precomputed[index]; } /*Elastic component (from Eq 17 in Adhikari et al, GMD 2015): */ if(computeelastic){ - int index=reCast(alpha/PI*(M-1)); - G_elastic[i] += G_elastic_precomputed[index]; - //values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i]; - Sgo[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i]; - } - } - - - /*Free ressources:*/ - if(computeelastic) xDelete(G_elastic); - if(computerigid) xDelete(G_rigid); + Sgo[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic_precomputed[index]; + } + } + return; Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24920) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24921) @@ -308,4 +308,5 @@ syn keyword cConstant SealevelriseEquatorialMoiEnum syn keyword cConstant SealevelriseFluidLoveEnum +syn keyword cConstant SealevelriseGRigidEnum syn keyword cConstant SealevelriseGElasticEnum syn keyword cConstant SealevelriseGeodeticEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24920) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24921) @@ -302,4 +302,5 @@ SealevelriseEquatorialMoiEnum, SealevelriseFluidLoveEnum, + SealevelriseGRigidEnum, SealevelriseGElasticEnum, SealevelriseGeodeticEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24920) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24921) @@ -310,4 +310,5 @@ case SealevelriseEquatorialMoiEnum : return "SealevelriseEquatorialMoi"; case SealevelriseFluidLoveEnum : return "SealevelriseFluidLove"; + case SealevelriseGRigidEnum : return "SealevelriseGRigid"; case SealevelriseGElasticEnum : return "SealevelriseGElastic"; case SealevelriseGeodeticEnum : return "SealevelriseGeodetic"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24920) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24921) @@ -316,4 +316,5 @@ else if (strcmp(name,"SealevelriseEquatorialMoi")==0) return SealevelriseEquatorialMoiEnum; else if (strcmp(name,"SealevelriseFluidLove")==0) return SealevelriseFluidLoveEnum; + else if (strcmp(name,"SealevelriseGRigid")==0) return SealevelriseGRigidEnum; else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum; else if (strcmp(name,"SealevelriseGeodetic")==0) return SealevelriseGeodeticEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum; else if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum; - else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum; + if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum; + else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum; else if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum; else if (strcmp(name,"SmbRunoffalti")==0) return SmbRunoffaltiEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum; else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum; - else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; + if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum; + else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum; else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; else if (strcmp(name,"Input")==0) return InputEnum; - else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum; + if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; + else if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum; else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum; else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum; - else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum; + if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum; + else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum; else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum; else if (strcmp(name,"SmbMeanULW")==0) return SmbMeanULWEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum; else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum; - else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum; + if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum; + else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum; else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum; else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum; @@ -997,9 +998,9 @@ else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum; - else if (strcmp(name,"BoolInput")==0) return BoolInputEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum; + if (strcmp(name,"BoolInput")==0) return BoolInputEnum; + else if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum; else if (strcmp(name,"IntInput2")==0) return IntInput2Enum; else if (strcmp(name,"BoolParam")==0) return BoolParamEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"IceVolumeAboveFloatationScaled")==0) return IceVolumeAboveFloatationScaledEnum; else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; - else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum; + if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum; + else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum; else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum; else if (strcmp(name,"Incremental")==0) return IncrementalEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"Profiler")==0) return ProfilerEnum; else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum; - else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum; + if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum; + else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum; else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum; else if (strcmp(name,"Regular")==0) return RegularEnum;