Changeset 21794

Timestamp:

07/13/17 17:02:28 (8 years ago)

Author:

seroussi

Message:

NEW: included improved treatment of melting under partially floating elements

File:

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

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

@@ -470 +470 @@
         /*Intermediaries */
         IssmDouble  Jdet,dt;
-        IssmDouble  ms,mb,gmb,fmb,thickness,phi;
+        IssmDouble  ms,mb,gmb,fmb,thickness;
+        IssmDouble  gllevelset,phi=1.;
         IssmDouble* xyz_list = NULL;
 
@@ -485 +486 @@
         Input* gmb_input           = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
         Input* fmb_input           = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
-        Input* groundedice_input   = element->GetInput(MaskGroundediceLevelsetEnum);              _assert_(groundedice_input);
+        Input* gllevelset_input   = element->GetInput(MaskGroundediceLevelsetEnum);              _assert_(gllevelset_input);
         Input* ms_input            = element->GetInput(SmbMassBalanceEnum);                       _assert_(ms_input);
         Input* thickness_input     = element->GetInput(ThicknessEnum);                            _assert_(thickness_input);
@@ -500 +501 @@
                 gmb_input->GetInputValue(&gmb,gauss);
                 fmb_input->GetInputValue(&fmb,gauss);
-                groundedice_input->GetInputValue(&phi,gauss);
+                gllevelset_input->GetInputValue(&phi,gauss);
                 thickness_input->GetInputValue(&thickness,gauss);
                 if(phi>0.) mb=gmb;
@@ -520 +521 @@
 
         /*Intermediaries */
-        IssmDouble  Jdet,dt;
-        IssmDouble  ms,mb,gmb,fmb,thickness,phi;
+        int         migration_style,point1;
+        bool        mainlyfloating;
+        IssmDouble  Jdet,fraction1,fraction2;
+        IssmDouble  dt,gllevelset;
+        IssmDouble  ms,mb,gmb,fmb,thickness,phi=1.;
         IssmDouble* xyz_list = NULL;
+        Gauss*      gauss     = NULL;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -537 +542 @@
         Input* fmb_input           = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
         Input* ms_input            = element->GetInput(SmbMassBalanceEnum);          _assert_(ms_input);
-        Input* groundedice_input   = element->GetInput(MaskGroundediceLevelsetEnum);             _assert_(groundedice_input);
+        Input* gllevelset_input   = element->GetInput(MaskGroundediceLevelsetEnum);             _assert_(gllevelset_input);
         Input* thickness_input     = element->GetInput(ThicknessEnum);                           _assert_(thickness_input);
 
+        /*Recover portion of element that is grounded*/
+        phi=element->GetGroundedPortion(xyz_list);
+        if(migration_style==SubelementMigration2Enum){
+                gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
+                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+           gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,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);
@@ -551 +566 @@
                 gmb_input->GetInputValue(&gmb,gauss);
                 fmb_input->GetInputValue(&fmb,gauss);
-                groundedice_input->GetInputValue(&phi,gauss);
+                gllevelset_input->GetInputValue(&gllevelset,gauss);
                 thickness_input->GetInputValue(&thickness,gauss);
-                if(phi>0) mb=gmb;
-                else mb=fmb;
+                if(migration_style==SubelementMigrationEnum){
+                        if (phi<0.00000001) mb=gmb;
+                        else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
+                }
+                else if(migration_style==SubelementMigration2Enum){
+                        if(gllevelset>0.) mb=gmb;
+                        else mb=fmb;
+                }
+                else if(migration_style==SubelementMigration3Enum){
+                        if (phi<0.00000001) mb=gmb;
+                        else mb=gmb;
+                }
+                else{
+                        if(gllevelset>0) mb=gmb;
+                        else mb=fmb;
+                }
 
                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];