Changeset 24944

Timestamp:

06/01/20 10:27:05 (5 years ago)

Author:

Eric.Larour

Message:

CHG: reorg and cleanup of the sea level cores, pushing as much as possible into the
geometry core.

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Tria.cpp (modified) (10 diffs)
• classes/FemModel.cpp (modified) (2 diffs)
• shared/Enum/Enum.vim (modified) (1 diff)
• shared/Enum/EnumDefinitions.h (modified) (1 diff)
• shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• shared/Enum/StringToEnumx.cpp (modified) (6 diffs)

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -5584 +5584 @@
         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<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
-        parameter->GetParameterValueByPointer(&dummypointer,&M);
-
+        parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M);
+
+        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
+        parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
+        
         /*allocate indices:*/
         indices=xNew<IssmDouble>(gsize);
+        G=xNewZeroInit<IssmDouble>(gsize);
+        
+        /*compute area:*/
+        area=GetAreaSpherical();
 
         /* Where is the centroid of this element?:{{{*/
@@ -5646 +5662 @@
         longe=longe/180*PI;
         /*}}}*/
-
+                
+        constant=3/rho_earth*area/eartharea;
+        
         for(int i=0;i<gsize;i++){
                 IssmDouble alpha;
@@ -5655 +5673 @@
                 delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
                 alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
-                indices[i]=reCast<int,IssmDouble>(alpha/PI*reCast<IssmDouble,int>(M-1));
+                indices[i]=alpha/PI*reCast<IssmDouble,int>(M-1);
+
+                /*Rigid earth gravitational perturbation: */
+                if(computerigid){
+                        G[i]+=G_rigid_precomputed[reCast<int,IssmDouble>(indices[i])];
+                }
+                if(computeelastic){
+                        G[i]+=G_elastic_precomputed[reCast<int,IssmDouble>(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 +5720 @@
 
         /*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 +5745 @@
         /*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 +5754 @@
         /*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 +5766 @@
     
         /*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<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
-        parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M);
-
-        if(computeelastic){
-                parameter = static_cast<DoubleVecParam*>(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 +5797 @@
         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 +5822 @@
                 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<gsize;i++){
-
-                        IssmDouble G_elastic=0;  //do not remove =0!
-
-                        /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
-                        /*lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
-                        delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
-                        alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
-                        index=reCast<int,IssmDouble>(alpha/PI*reCast<IssmDouble,int>(M-1));*/
-                        index=reCast<int,IssmDouble>(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<gsize;i++) Sgi[i]+=G[i]*I;
 
         /*Assign output pointer:*/
         _assert_(!xIsNan<IssmDouble>(eustatic));
-        _assert_(!xIsInf<IssmDouble>(eustatic));
         *peustatic=eustatic;
         return;
@@ -5998 +5994 @@
         /*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;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[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<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
-        parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M);
-
-        if(computeelastic){
-                /*recover elastic green function:*/
-                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
-                parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
-        }
-
-        
-        for(int i=0;i<gsize;i++){
-
-                /*Compute alpha angle between centroid and current vertex : */
-                /*lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
-                delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
-                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
-                index=reCast<int,IssmDouble>(alpha/PI*(M-1));*/
-                        
-                index=reCast<int,IssmDouble>(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;i<gsize;i++) Sgo[i]+=G[i]*S;
 
         return;

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -4745 +4745 @@
         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:*/
@@ -4754 +4761 @@
         RSLg_old=pRSLg_old->ToMPISerial();
 
-        /*Call the sea level rise core: */
-        for(int i=0;i<elements->Size();i++){
-                Element* element=xDynamicCast<Element*>(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;i<elements->Size();i++){
+                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                        element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks, latitude,longitude,radius);
+                }
         }
         pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL);

issm/trunk-jpl/src/c/shared/Enum/Enum.vim

@@ -685 +685 @@
 syn keyword cConstant SealevelriseHydroRateEnum
 syn keyword cConstant SealevelriseIndicesEnum
+syn keyword cConstant SealevelriseGEnum
 syn keyword cConstant SedimentHeadEnum
 syn keyword cConstant SedimentHeadOldEnum

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -682 +682 @@
         SealevelriseHydroRateEnum,
         SealevelriseIndicesEnum,
+        SealevelriseGEnum,
         SedimentHeadEnum,
         SedimentHeadOldEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -687 +687 @@
                 case SealevelriseHydroRateEnum : return "SealevelriseHydroRate";
                 case SealevelriseIndicesEnum : return "SealevelriseIndices";
+                case SealevelriseGEnum : return "SealevelriseG";
                 case SedimentHeadEnum : return "SedimentHead";
                 case SedimentHeadOldEnum : return "SedimentHeadOld";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -702 +702 @@
               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 +752 @@
               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 +875 @@
               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 +998 @@
               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 +1121 @@
               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 +1244 @@
               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;