Changeset 24117

Timestamp:

08/05/19 09:38:22 (6 years ago)

Author:

tsantos

Message:

CHG: extending subelement melt parameterizations to discontinuous Galerkin formulation

File:

• issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (modified) (4 diffs)

issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

@@ -694 +694 @@
 
         /*Intermediaries */
+        int melt_style, point1;
+        bool mainlyfloating;
+        IssmDouble  fraction1,fraction2,gllevelset;
         IssmDouble  Jdet,dt;
         IssmDouble  ms,mb,gmb,fmb,thickness,phi=1.;
@@ -708 +711 @@
         element->GetVerticesCoordinates(&xyz_list);
         element->FindParam(&dt,TimesteppingTimeStepEnum);
+        element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum);
         Input* gmb_input        = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
         Input* fmb_input        = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
@@ -714 +718 @@
         Input* thickness_input  = element->GetInput(ThicknessEnum);                           _assert_(thickness_input);
 
+   /*Recover portion of element that is grounded*/
+   Gauss* gauss=NULL;
+   phi=element->GetGroundedPortion(xyz_list);
+   if(melt_style==SubelementMelt2Enum){
+      element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+      gauss = element->NewGauss(point1,fraction1,fraction2,3);
+   }
+   else{
+      gauss = element->NewGauss(3);
+   }
+
         /* Start  looping on the number of gaussian points: */
-        Gauss* gauss=element->NewGauss(2);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
                 gauss->GaussPoint(ig);
@@ -725 +739 @@
                 gmb_input->GetInputValue(&gmb,gauss);
                 fmb_input->GetInputValue(&fmb,gauss);
-                gllevelset_input->GetInputValue(&phi,gauss);
+                gllevelset_input->GetInputValue(&gllevelset,gauss);
                 thickness_input->GetInputValue(&thickness,gauss);
 
-                if(phi>0.) mb=gmb;
-                else mb=fmb;
+                if(melt_style==SubelementMelt1Enum){
+         if (phi>0.999999999) mb=gmb;
+         else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
+      }
+      else if(melt_style==SubelementMelt2Enum){
+         if(gllevelset>0.) mb=gmb;
+         else mb=fmb;
+      }
+      else if(melt_style==NoMeltOnPartiallyFloatingEnum){
+         if (phi<0.00000001) mb=fmb;
+         else mb=gmb;
+      }
+      else if(melt_style==FullMeltOnPartiallyFloatingEnum){
+         if (phi<0.99999999) mb=fmb;
+         else mb=gmb;
+      }
+      else  _error_("melt interpolation "<<EnumToStringx(melt_style)<<" not implemented yet");
 
                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];