Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24943) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24944) @@ -5584,27 +5584,43 @@ bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble area; + IssmDouble area,eartharea; IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) - IssmDouble rho; + IssmDouble rho_earth; IssmDouble late,longe,re; IssmDouble lati,longi,ri; + IssmDouble constant; + + IssmDouble* G=NULL; + IssmDouble* G_elastic_precomputed=NULL; + IssmDouble* G_rigid_precomputed=NULL; /*elastic green function:*/ IssmDouble* indices=NULL; int M; - IssmDouble* dummypointer=NULL; /*Computational flags:*/ bool computerigid = true; - - /*how many dofs are we working with here? */ + bool computeelastic = true; + + /*recover parameters: */ + rho_earth=FindParam(MaterialsEarthDensityEnum); + this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); + this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); - - /*recover size of indexed tables:*/ + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); + + /*recover precomputed green function kernels:*/ DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer(&dummypointer,&M); - + parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M); + + parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M); + /*allocate indices:*/ indices=xNew(gsize); + G=xNewZeroInit(gsize); + + /*compute area:*/ + area=GetAreaSpherical(); /* Where is the centroid of this element?:{{{*/ @@ -5646,5 +5662,7 @@ longe=longe/180*PI; /*}}}*/ - + + constant=3/rho_earth*area/eartharea; + for(int i=0;i(alpha/PI*reCast(M-1)); + indices[i]=alpha/PI*reCast(M-1); + + /*Rigid earth gravitational perturbation: */ + if(computerigid){ + G[i]+=G_rigid_precomputed[reCast(indices[i])]; + } + if(computeelastic){ + G[i]+=G_elastic_precomputed[reCast(indices[i])]; + } + G[i]=G[i]*constant; } /*Add in inputs:*/ this->inputs2->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); - - /*Another quantity, area:*/ - this->inputs2->SetDoubleInput(AreaEnum,this->lid,GetAreaSpherical()); + this->inputs2->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); + this->inputs2->SetDoubleInput(AreaEnum,this->lid,area); + + /*Free allocations:*/ + xDelete(indices); + xDelete(G); return; @@ -5690,20 +5720,12 @@ /*diverse:*/ - int gsize,dummy; - bool spherical=true; - IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble area,eartharea; + int gsize; + IssmDouble area; + IssmDouble phi=1.0; //WARNING: do not touch this, default is entire elemnt contributes eustatic IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) - IssmDouble rho; - IssmDouble late,longe,re; - IssmDouble lati,longi,ri; bool notfullygrounded=false; /*elastic green function:*/ - int index; - IssmDouble* indices=NULL; - IssmDouble* G_elastic_precomputed=NULL; - IssmDouble* G_rigid_precomputed=NULL; - int M; + IssmDouble* G=NULL; /*ice properties: */ @@ -5723,5 +5745,7 @@ /*early return if we are not on an ice cap:*/ if(!masks->isiceonly[this->lid]){ + #ifdef _ISSM_DEBUG_ constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); + #endif *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! return; @@ -5730,5 +5754,8 @@ /*early return if we are fully floating:*/ if (masks->isfullyfloating[this->lid]){ - constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); + constant=0; + #ifdef _ISSM_DEBUG_ + this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); + #endif *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! return; @@ -5739,45 +5766,30 @@ /*Inform mask: */ - constant=1; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); + constant=1; + #ifdef _ISSM_DEBUG_ + this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); + #endif /*recover material parameters: */ rho_ice=FindParam(MaterialsRhoIceEnum); rho_water=FindParam(MaterialsRhoSeawaterEnum); - rho_earth=FindParam(MaterialsEarthDensityEnum); /*recover love numbers and computational flags: */ - this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); - this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum); - /*recover earth area: */ - this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); - - /*recover precomputed green function kernels:*/ - DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M); - - if(computeelastic){ - parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M); - } - - /*how many dofs are we working with here? */ - this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); - - /*retrieve indices:*/ - if(computerigid){this->inputs2->GetArrayPtr(SealevelriseIndicesEnum,this->lid,&indices,&dummy); _assert_(dummy==gsize);} + /*retrieve precomputed G:*/ + this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ this->GetInput2Value(&area,AreaEnum); - if(notfullygrounded){ - IssmDouble phi=0; + /*factor to reduce area if we are not fully grounded:*/ + if(notfullygrounded){ /*{{{*/ IssmDouble xyz_list[NUMVERTICES][3]; ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias - area*=phi; //scale our load by the fraction of grounded area. - } + } + /*}}}*/ /*Compute ice thickness: */ @@ -5785,5 +5797,5 @@ if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); - /*If we are fully grounded, take the average over the element: */ + /*If we are fully grounded, take the average over the element: {{{*/ if(!notfullygrounded)deltathickness_input->GetInputAverage(&I); else{ @@ -5810,35 +5822,19 @@ delete gauss; } + /*}}}*/ /*Compute eustatic compoent:*/ _assert_(oceanarea>0.); if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 - eustatic += rho_ice*area*I/(oceanarea*rho_water); - - if(computeelastic | computerigid){ - IssmDouble alpha; - IssmDouble delPhi,delLambda; - for(int i=0;i(alpha/PI*reCast(M-1));*/ - index=reCast(indices[i]); - - //Elastic component (from Eq 17 in Adhikari et al, GMD 2015) - 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_precomputed[index]+G_elastic); - } - } + eustatic += rho_ice*area*phi*I/(oceanarea*rho_water); + + /*convert from m to kg/m^2:*/ + I=I*rho_ice*phi; + + /*convolve:*/ + for(int i=0;i(eustatic)); - _assert_(!xIsInf(eustatic)); *peustatic=eustatic; return; @@ -5998,98 +5994,40 @@ /*diverse:*/ int gsize,dummy; - bool spherical=true; - IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble area,eartharea; IssmDouble S; //change in water water level(Farrel and Clarke, Equ. 4) - IssmDouble late,longe; - IssmDouble lati,longi,ri; - IssmDouble minlong=400; - IssmDouble maxlong=-20; IssmDouble constant=0; - - /*precomputed elastic green functions:*/ - IssmDouble* G_rigid_precomputed = NULL; - IssmDouble* G_elastic_precomputed = NULL; - int M; - int index; - IssmDouble* indices=NULL; - IssmDouble alpha; - IssmDouble delPhi,delLambda; + IssmDouble rho_water; + IssmDouble* G=NULL; /*optimization:*/ bool store_green_functions=false; - - /*ice properties: */ - IssmDouble rho_ice,rho_water,rho_earth; - - /*Computational flags:*/ - bool computerigid = true; - bool computeelastic= true; + + rho_water=FindParam(MaterialsRhoSeawaterEnum); /*early return if we are not on the ocean:*/ if (!masks->isoceanin[this->lid]){ - constant=0; this->AddInput2(SealevelEustaticOceanMaskEnum,&constant,P0Enum); + constant=0; + #ifdef _ISSM_DEBUG_ + this->AddInput2(SealevelEustaticOceanMaskEnum,&constant,P0Enum); + #endif return; } - constant=1; this->AddInput2(SealevelEustaticOceanMaskEnum,&constant,P0Enum); - - /*recover computational flags: */ - this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); - this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); - - /*recover earth area: */ - this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); - - /*early return if rigid or elastic not requested:*/ - if(!computerigid && !computeelastic) return; - - /*recover material parameters: */ - rho_water=FindParam(MaterialsRhoSeawaterEnum); - rho_earth=FindParam(MaterialsEarthDensityEnum); + constant=1; + #ifdef _ISSM_DEBUG_ + this->AddInput2(SealevelEustaticOceanMaskEnum,&constant,P0Enum); + #endif /*how many dofs are we working with here? */ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); - - /*retrieve indices:*/ - if(computerigid){this->inputs2->GetArrayPtr(SealevelriseIndicesEnum,this->lid,&indices,&dummy); _assert_(dummy==gsize);} + + /*retrieve precomputed G:*/ + this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&dummy); _assert_(dummy==gsize); /*From Sg_old, recover water sea level rise:*/ S=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; - /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ - this->GetInput2Value(&area,AreaEnum); - - /*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:*/ - parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M); - } - - - for(int i=0;i(alpha/PI*(M-1));*/ - - index=reCast(indices[i]); - - /*Rigid earth gravitational perturbation: */ - if(computerigid){ - 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){ - Sgo[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic_precomputed[index]; - } - } - + /*convert to SI: */ + S=S*rho_water; + + for(int i=0;iparameters->FindParam(&computerigid,SealevelriseRigidEnum); + this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); + /*Initialize temporary vector that will be used to sum eustatic components on all local elements, prior * to assembly:*/ @@ -4754,8 +4761,10 @@ RSLg_old=pRSLg_old->ToMPISerial(); - /*Call the sea level rise core: */ - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks, latitude,longitude,radius); + /*Call the sea level rise non-eustatic core only if required: */ + if(computerigid | computeelastic){ + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks, latitude,longitude,radius); + } } pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL); Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24943) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24944) @@ -685,4 +685,5 @@ syn keyword cConstant SealevelriseHydroRateEnum syn keyword cConstant SealevelriseIndicesEnum +syn keyword cConstant SealevelriseGEnum syn keyword cConstant SedimentHeadEnum syn keyword cConstant SedimentHeadOldEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24943) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24944) @@ -682,4 +682,5 @@ SealevelriseHydroRateEnum, SealevelriseIndicesEnum, + SealevelriseGEnum, SedimentHeadEnum, SedimentHeadOldEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24943) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24944) @@ -687,4 +687,5 @@ case SealevelriseHydroRateEnum : return "SealevelriseHydroRate"; case SealevelriseIndicesEnum : return "SealevelriseIndices"; + case SealevelriseGEnum : return "SealevelriseG"; case SedimentHeadEnum : return "SedimentHead"; case SedimentHeadOldEnum : return "SedimentHeadOld"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24943) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24944) @@ -702,4 +702,5 @@ else if (strcmp(name,"SealevelriseHydroRate")==0) return SealevelriseHydroRateEnum; else if (strcmp(name,"SealevelriseIndices")==0) return SealevelriseIndicesEnum; + else if (strcmp(name,"SealevelriseG")==0) return SealevelriseGEnum; else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum; else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; - else if (strcmp(name,"SmbMassBalanceClimate")==0) return SmbMassBalanceClimateEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; + if (strcmp(name,"SmbMassBalanceClimate")==0) return SmbMassBalanceClimateEnum; + else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum; else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"Outputdefinition10")==0) return Outputdefinition10Enum; else if (strcmp(name,"Outputdefinition11")==0) return Outputdefinition11Enum; - else if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum; + if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum; + else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum; else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum; else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum; @@ -997,9 +998,9 @@ else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum; else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum; - else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum; + if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum; + else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum; else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"Hydrologyshakti")==0) return HydrologyshaktiEnum; else if (strcmp(name,"Hydrologyshreve")==0) return HydrologyshreveEnum; - else if (strcmp(name,"IceMass")==0) return IceMassEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"IceMassScaled")==0) return IceMassScaledEnum; + if (strcmp(name,"IceMass")==0) return IceMassEnum; + else if (strcmp(name,"IceMassScaled")==0) return IceMassScaledEnum; else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum; else if (strcmp(name,"IceVolumeAboveFloatationScaled")==0) return IceVolumeAboveFloatationScaledEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"Paterson")==0) return PatersonEnum; else if (strcmp(name,"Pengrid")==0) return PengridEnum; - else if (strcmp(name,"Penpair")==0) return PenpairEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"Penta")==0) return PentaEnum; + if (strcmp(name,"Penpair")==0) return PenpairEnum; + else if (strcmp(name,"Penta")==0) return PentaEnum; else if (strcmp(name,"PentaInput")==0) return PentaInputEnum; else if (strcmp(name,"Profiler")==0) return ProfilerEnum;