Changeset 24968

Timestamp:

06/04/20 19:56:33 (5 years ago)

Author:

Eric.Larour

Message:

CHG: bottom pressure fingerprints computations.

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Tria.cpp (modified) (8 diffs)
• Tria.h (modified) (1 diff)

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

@@ -5751 +5751 @@
 void    Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble oceanarea){ /*{{{*/
 
-        /*Computational flags:*/
         int bp_compute_fingerprints= 0;
-
-        /*some paramters first: */
+                
+        /*Compute bottom pressure contribution from ocean if requested:*/
         this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
-
-        if(!masks->isoceanin[this->lid]){
-                /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested
-                 *bottom pressure fingerprints:*/
-                if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,oceanarea);
-        }
-        //if(!IsIceInElement()){
-                /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/
-                this->SealevelriseEustaticIce(Sgi,peustatic,masks, oceanarea);
-        //}
+        if(bp_compute_fingerprints)this->SealevelriseBottomPressure(Sgi,masks);
+
+        /*Compute eustatic ice contribution to sea level rise: */
+        this->SealevelriseEustaticIce(Sgi,peustatic,masks, oceanarea);
 
 }
@@ -5777 +5770 @@
         IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
         bool notfullygrounded=false;
+        bool scaleoceanarea= false;
 
         /*elastic green function:*/
@@ -5782 +5776 @@
 
         /*ice properties: */
-        IssmDouble rho_ice,rho_water,rho_earth;
+        IssmDouble rho_ice,rho_water;
 
         /*constants:*/
@@ -5789 +5783 @@
         /*Initialize eustatic component: do not skip this step :):*/
         IssmDouble eustatic = 0.;
-
-        /*Computational flags:*/
-        bool computerigid = true;
-        bool computeelastic= true;
-        bool scaleoceanarea= false;
 
         /*early return if we are not on an ice cap:*/
@@ -5827 +5816 @@
         rho_water=FindParam(MaterialsRhoSeawaterEnum);
 
-        /*recover love numbers and computational flags: */
+        /*recover ocean area scaling: */
         this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum);
 
@@ -5836 +5825 @@
         this->GetInput2Value(&area,AreaEnum);
 
-        /*factor to reduce area if we are not fully grounded:*/
-        if(notfullygrounded){ /*{{{*/
+        /*Compute fraction of the element that is 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
-        }
-        /*}}}*/
-
-        /*Compute ice thickness: */
+        } 
+        else phi=1.0;
+
+        /*Retrieve ice thickness at vertices: */
         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: {{{*/
+        /*/Average ice thickness over grounded area of the element only: {{{*/
         if(!notfullygrounded)deltathickness_input->GetInputAverage(&I);
         else{
@@ -5876 +5865 @@
         /*}}}*/
 
-        /*Compute eustatic compoent:*/
+        /*Compute eustatic component:*/
         _assert_(oceanarea>0.);
         if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
@@ -5893 +5882 @@
 }
 /*}}}*/
-void    Tria::SealevelriseEustaticBottomPressure(IssmDouble* Sgi,IssmDouble* peustatic, IssmDouble oceanarea){ /*{{{*/
+void    Tria::SealevelriseBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/
 
         /*diverse:*/
         int gsize;
-        bool spherical=true;
-        IssmDouble llr_list[NUMVERTICES][3];
-        IssmDouble area,planetarea;
-        IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
-        IssmDouble rho;
-        IssmDouble late,longe,re;
-        IssmDouble lati,longi,ri;
+        IssmDouble area;
+        IssmDouble BP;  //change in bottom pressure (Farrel and Clarke, Equ. 4)
 
         /*elastic green function:*/
-        IssmDouble* G_elastic_precomputed=NULL;
-        int         M;
-
-        /*ice properties: */
-        IssmDouble rho_water,rho_earth;
-
-        /*constants:*/
-        IssmDouble constant=0;
-
-        /*Initialize eustatic component: do not skip this step :):*/
-        IssmDouble eustatic = 0.;
-
-        /*Computational flags:*/
-        bool computerigid = true;
-        bool computeelastic= true;
-        bool scaleoceanarea= false;
-        bool bp_compute_fingerprints= false;
-
-        /*we are here to compute fingerprints originating fromn bottom pressure loads:*/
+        IssmDouble* G=NULL;
+        
+        /*water properties: */
+        IssmDouble rho_water;
+
+        /*we are here to compute fingerprints originating fromn bottom pressure changes:*/
+        if(!masks->isoceanin[this->lid])return;
 
         /*Inform mask: */
+        #ifdef _ISSM_DEBUG_
         constant=1; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
+        #endif
 
         /*recover material parameters: */
         rho_water=FindParam(MaterialsRhoSeawaterEnum);
-        rho_earth=FindParam(MaterialsEarthDensityEnum);
-
-        /*recover earth area: */
-        this->parameters->FindParam(&planetarea,SealevelPlanetAreaEnum);
-
-        /*recover love numbers and computational flags: */
-        this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
-        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
-        this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum);
-
-        /*recover elastic green function:*/
-        if(computeelastic){
-                DoubleVecParam* 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);
-
-        /* Where is the centroid of this element?:{{{*/
-
-        /*retrieve coordinates: */
-        ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
-
-        IssmDouble minlong=400;
-        IssmDouble maxlong=-20;
-        for (int i=0;i<NUMVERTICES;i++){
-                llr_list[i][0]=(90-llr_list[i][0]);
-                if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
-                if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
-                if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
-        }
-        if(minlong==0 && maxlong>180){
-                if (llr_list[0][1]==0)llr_list[0][1]=360;
-                if (llr_list[1][1]==0)llr_list[1][1]=360;
-                if (llr_list[2][1]==0)llr_list[2][1]=360;
-        }
-
-        // correction at the north pole
-        if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
-        if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
-        if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
-
-        //correction at the south pole
-        if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
-        if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
-        if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
-
-        late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
-        longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
-
-        late=90-late;
-        if(longe>180)longe=(longe-180)-180;
-
-        late=late/180*PI;
-        longe=longe/180*PI;
-        /*}}}*/
+
+        /*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);
 
-        /*Compute bottom pressure change: */
+        /*Retrieve bottom pressure change and average over the element: */
         Input2* bottompressure_change_input=this->GetInput2(DslSeaWaterPressureChangeAtSeaFloor);
         if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level rise fingerprint!");
-
-        /*If we are fully grounded, take the average over the element: */
-        bottompressure_change_input->GetInputAverage(&I);
-
-        /*Compute eustatic compoent:*/
-        _assert_(oceanarea>0.);
-        if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
-
-        /*We do not need to add the bottom pressure component to the eustatic value: */
-        eustatic += 0;
-        IssmDouble* latitude=NULL; //NOT GOING TO WORK!
-        IssmDouble* longitude=NULL; //NOT GOING TO WORK!
-
-        if(computeelastic | computerigid){
-                IssmDouble alpha;
-                IssmDouble delPhi,delLambda;
-                for(int i=0;i<gsize;i++){
-
-                        IssmDouble G_rigid=0;  //do not remove =0!
-                        IssmDouble G_elastic=0;  //do not remove =0!
-
-                        /*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)));
-
-                        //Rigid earth gravitational perturbation:
-                        if(computerigid)G_rigid=1.0/2.0/sin(alpha/2.0);
-
-                        //Elastic component  (from Eq 17 in Adhikari et al, GMD 2015)
-                        if(computeelastic){
-                                int index=reCast<int,IssmDouble>(alpha/PI*reCast<IssmDouble,int>(M-1));
-                                G_elastic += G_elastic_precomputed[index];
-                        }
-
-                        /*Add all components to the pSgi or pSgo solution vectors:*/
-                        Sgi[i]+=3*rho_water/rho_earth*area/planetarea*I*(G_rigid+G_elastic);
-                }
-        }
-
-        /*Assign output pointer:*/
-        _assert_(!xIsNan<IssmDouble>(eustatic));
-        _assert_(!xIsInf<IssmDouble>(eustatic));
-        *peustatic=eustatic;
+        bottompressure_change_input->GetInputAverage(&BP);
+
+        /*convert from m to kg/m^2:*/
+        BP=BP*rho_water;
+
+        /*convolve:*/
+        for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*BP;
+
         return;
 }

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

@@ -169 +169 @@
                 void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble oceanarea);
                 void    SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble oceanarea);
-                void    SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble oceanarea);
+                void    SealevelriseBottomPressure(IssmDouble* Sgi,SealevelMasks* masks);
                 void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks);
                 void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks);