Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 24943) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 24944) @@ -4744,6 +4744,13 @@ int* indices = NULL; int gsize; + bool computerigid = true; + bool computeelastic= true; + + /*recover computational flags: */ + this->parameters->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:*/ gsize = this->nodes->NumberOfDofs(GsetEnum); @@ -4753,10 +4760,12 @@ /*Serialize vectors from previous iteration:*/ 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); pRSLgo->Assemble(); Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24943) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24944) @@ -5583,29 +5583,45 @@ int gsize; 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; + bool computeelastic = true; - /*how many dofs are we working with here? */ + /*recover parameters: */ + rho_earth=FindParam(MaterialsEarthDensityEnum); + this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); + this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); - /*recover size of indexed tables:*/ + /*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?:{{{*/ @@ -5645,7 +5661,9 @@ late=late/180*PI; 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); + this->inputs2->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); + this->inputs2->SetDoubleInput(AreaEnum,this->lid,area); - /*Another quantity, area:*/ - this->inputs2->SetDoubleInput(AreaEnum,this->lid,GetAreaSpherical()); + /*Free allocations:*/ + xDelete(indices); + xDelete(G); return; } @@ -5689,22 +5719,14 @@ void Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){ /*{{{*/ /*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: */ IssmDouble rho_ice,rho_water,rho_earth; @@ -5722,7 +5744,9 @@ /*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; } @@ -5729,7 +5753,10 @@ /*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; } @@ -5738,53 +5765,38 @@ if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on. /*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); + /*retrieve precomputed G:*/ + this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); - /*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);} - /*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: */ Input2* deltathickness_input=this->GetInput2(SealevelriseDeltathicknessEnum); 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{ IssmDouble total_weight=0; @@ -5809,37 +5821,21 @@ I=I/total_weight; 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); + eustatic += rho_ice*area*phi*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); - } - } - /*Assign output pointer:*/ _assert_(!xIsNan(eustatic)); - _assert_(!xIsInf(eustatic)); *peustatic=eustatic; return; } @@ -5997,101 +5993,43 @@ /*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; + IssmDouble rho_water; + IssmDouble* G=NULL; - /*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; - /*optimization:*/ bool store_green_functions=false; + + rho_water=FindParam(MaterialsRhoSeawaterEnum); - /*ice properties: */ - IssmDouble rho_ice,rho_water,rho_earth; - - /*Computational flags:*/ - bool computerigid = true; - bool computeelastic= true; - /*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); + constant=1; + #ifdef _ISSM_DEBUG_ + this->AddInput2(SealevelEustaticOceanMaskEnum,&constant,P0Enum); + #endif - /*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); - /*how many dofs are we working with here? */ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + + /*retrieve precomputed G:*/ + this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&dummy); _assert_(dummy==gsize); - /*retrieve indices:*/ - if(computerigid){this->inputs2->GetArrayPtr(SealevelriseIndicesEnum,this->lid,&indices,&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); + /*convert to SI: */ + S=S*rho_water; - /*recover rigid and elastic green functions:*/ - DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M); + for(int i=0;i(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]; - } - } - - return; } /*}}}*/ Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24943) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24944) @@ -681,6 +681,7 @@ SealevelriseSpcthicknessEnum, SealevelriseHydroRateEnum, SealevelriseIndicesEnum, + SealevelriseGEnum, SedimentHeadEnum, SedimentHeadOldEnum, SedimentHeadSubstepEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24943) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24944) @@ -686,6 +686,7 @@ case SealevelriseSpcthicknessEnum : return "SealevelriseSpcthickness"; case SealevelriseHydroRateEnum : return "SealevelriseHydroRate"; case SealevelriseIndicesEnum : return "SealevelriseIndices"; + case SealevelriseGEnum : return "SealevelriseG"; case SedimentHeadEnum : return "SedimentHead"; case SedimentHeadOldEnum : return "SedimentHeadOld"; case SedimentHeadSubstepEnum : return "SedimentHeadSubstep"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24943) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24944) @@ -701,6 +701,7 @@ else if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum; 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; else if (strcmp(name,"SedimentHeadSubstep")==0) return SedimentHeadSubstepEnum; @@ -750,11 +751,11 @@ else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum; 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; else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum; @@ -873,11 +874,11 @@ else if (strcmp(name,"Outputdefinition1")==0) return Outputdefinition1Enum; 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; else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum; @@ -996,11 +997,11 @@ else if (strcmp(name,"BalancevelocityAnalysis")==0) return BalancevelocityAnalysisEnum; 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; else if (strcmp(name,"BoolInput")==0) return BoolInputEnum; @@ -1119,11 +1120,11 @@ else if (strcmp(name,"Hydrologypism")==0) return HydrologypismEnum; 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; else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; @@ -1242,11 +1243,11 @@ else if (strcmp(name,"P2xP4")==0) return P2xP4Enum; 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; else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum; Index: ../trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24943) +++ ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24944) @@ -684,6 +684,7 @@ syn keyword cConstant SealevelriseSpcthicknessEnum syn keyword cConstant SealevelriseHydroRateEnum syn keyword cConstant SealevelriseIndicesEnum +syn keyword cConstant SealevelriseGEnum syn keyword cConstant SedimentHeadEnum syn keyword cConstant SedimentHeadOldEnum syn keyword cConstant SedimentHeadSubstepEnum