Changeset 26158

Timestamp:

03/29/21 09:14:55 (4 years ago)

Author:

Eric.Larour

Message:

CHG: fixed old code issues with masks in sea level core.

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (4 diffs)

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

@@ -5534 +5534 @@
                 dI_list[2] = +4*M_PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
         }
-        else if(masks->isiceonly[this->lid]){
-                IssmDouble rho_ice, I;
-
-                /*recover material parameters: */
+        if(masks->isiceonly[this->lid]){
+                
+                IssmDouble rho_ice,I;
+                bool computerigid= false;
+        
+                bool notfullygrounded=false;
+        
+                bool scaleoceanarea= false;
+                int  glfraction=1;
+                IssmDouble phi=1.0;
+                /*{{{*/
+        
+                
+
+                if (masks->isfullyfloating[this->lid]){
+                        I=0;
+                        this->AddInput(EffectivePressureEnum,&I,P0Enum);
+                        return;
+                }
+                        
+
+        
+                if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on.
+
+                /*recover some parameters:*/
                 rho_ice=FindParam(MaterialsRhoIceEnum);
-
-                /*Compute ice thickness change: */
+                this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
+                this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
+                this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum);
+
+                /*Get area of element: precomputed in the sealevelchange_geometry:*/
+                this->Element::GetInputValue(&area,AreaEnum);
+
+                /*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
+                        if(glfraction==0)phi=1;
+                        #ifdef _ISSM_DEBUG_
+                        this->AddInput(SealevelBarystaticMaskEnum,&phi,P0Enum);
+                        #endif
+                }
+                else phi=1.0;
+
+                /*Retrieve surface load for ice: */
                 Input* deltathickness_input=this->GetInput(DeltaIceThicknessEnum);
-                if (!deltathickness_input)_error_("DeltaIceThicknessEnum delta ice thickness input needed to compute sea level change!");
-                deltathickness_input->GetInputAverage(&I);
+                if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!");
+
+                /*/Average ice thickness over grounded area of the element only: {{{*/
+                if(!notfullygrounded)deltathickness_input->GetInputAverage(&I);
+                else{
+                        IssmDouble total_weight=0;
+                        bool mainlyfloating = true;
+                        int         point1;
+                        IssmDouble  fraction1,fraction2;
+
+                        /*Recover portion of element that is grounded*/
+                        this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+                        Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
+
+                        /* Start  looping on the number of gaussian points and average over these gaussian points: */
+                        total_weight=0;
+                        I=0;
+                        while(gauss->next()){
+                                IssmDouble Ig=0;
+                                deltathickness_input->GetInputValue(&Ig,gauss);
+                                I+=Ig*gauss->weight;
+                                total_weight+=gauss->weight;
+                        }
+                        I=I/total_weight;
+                        delete gauss;
+                }
+                /*}}}*/
+
+                /*}}}*/
+                /*convert to kg/m^2*/
+                I=I*phi;
 
                 dI_list[0] = -4*M_PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
@@ -6053 +6122 @@
         /*diverse:*/
         int gsize;
-        IssmDouble I, S, BP;            //change in relative ice thickness and sea level
+        IssmDouble I=0;
+        IssmDouble S=0;
+        IssmDouble BP=0;
         IssmDouble rho_ice,rho_water;
         int horiz;
@@ -6065 +6136 @@
         /*computational flags:*/
         bool computeelastic= false;
+        bool computerigid= false;
+        bool notfullygrounded=false;
+        bool scaleoceanarea= false;
+        int  glfraction=1;
+        IssmDouble area;
+        IssmDouble phi=1.0;
 
         /*early return if we are not on the ocean or on an ice cap:*/
@@ -6116 +6193 @@
                 }
         }
-        else if (masks->isiceonly[this->lid]){
-
-                /*Compute ice thickness change: */
+        if (masks->isiceonly[this->lid]){
+
+                if (masks->isfullyfloating[this->lid]) return;
+
+        
+                if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on.
+
+                /*recover some parameters:*/
+                rho_ice=FindParam(MaterialsRhoIceEnum);
+                rho_water=FindParam(MaterialsRhoSeawaterEnum);
+                this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
+                this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
+                this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum);
+
+                /*Get area of element: precomputed in the sealevelchange_geometry:*/
+                this->Element::GetInputValue(&area,AreaEnum);
+
+                /*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
+                        if(glfraction==0)phi=1;
+                        #ifdef _ISSM_DEBUG_
+                        this->AddInput(SealevelBarystaticMaskEnum,&phi,P0Enum);
+                        #endif
+                }
+                else phi=1.0;
+
+                /*Retrieve surface load for ice: */
                 Input* deltathickness_input=this->GetInput(DeltaIceThicknessEnum);
-                if (!deltathickness_input)_error_("DeltaIceThicknessEnum delta thickness input needed to compute sea level change!");
-                deltathickness_input->GetInputAverage(&I);
+                if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!");
+
+                /*/Average ice thickness over grounded area of the element only: {{{*/
+                if(!notfullygrounded)deltathickness_input->GetInputAverage(&I);
+                else{
+                        IssmDouble total_weight=0;
+                        bool mainlyfloating = true;
+                        int         point1;
+                        IssmDouble  fraction1,fraction2;
+
+                        /*Recover portion of element that is grounded*/
+                        this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+                        Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
+
+                        /* Start  looping on the number of gaussian points and average over these gaussian points: */
+                        total_weight=0;
+                        I=0;
+                        while(gauss->next()){
+                                IssmDouble Ig=0;
+                                deltathickness_input->GetInputValue(&Ig,gauss);
+                                I+=Ig*gauss->weight;
+                                total_weight+=gauss->weight;
+                        }
+                        I=I/total_weight;
+                        delete gauss;
+                }
+                /*}}}*/
 
                 /*convert to kg/m^2*/
-                I=I*rho_ice;
-
+                I=I*rho_ice*phi;
+                
                 for(int i=0;i<gsize;i++){
                         Up[i]+=I*GU[i];