Changeset 11355

Timestamp:

02/07/12 15:49:52 (13 years ago)

Author:

seroussi

Message:

got confused between thermal and enthalpy model ...

File:

• issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (modified) (7 diffs)

issm/trunk-jpl/src/c/objects/Elements/Penta.cpp

@@ -3854 +3854 @@
         double     rho_ice,heatcapacity,geothermalflux_value;
         double     basalfriction,alpha2,vx,vy;
-        double     scalar;
+        double     scalar,enthalpy,enthalpyup;
+        double     pressure,pressureup;
         double     basis[NUMVERTICES];
         Friction*  friction=NULL;
         GaussPenta* gauss=NULL;
+        GaussPenta* gaussup=NULL;
 
         /* Geothermal flux on ice sheet base and basal friction */
@@ -3875 +3877 @@
         Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
         Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
+        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
+        Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
         Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
 
@@ -3882 +3886 @@
         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
         gauss=new GaussPenta(0,1,2,2);
+        gaussup=new GaussPenta(3,4,5,2);
         for(ig=gauss->begin();ig<gauss->end();ig++){
 
                 gauss->GaussPoint(ig);
+                gaussup->GaussPoint(ig);
 
                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
                 GetNodalFunctionsP1(&basis[0], gauss);
 
-                geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
-                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
-                vx_input->GetInputValue(&vx,gauss);
-                vy_input->GetInputValue(&vy,gauss);
-                basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
-                
-                scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
-                if(dt) scalar=dt*scalar;
-
-                for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
+                enthalpy_input->GetInputValue(&enthalpy,gauss);
+                pressure_input->GetInputValue(&pressure,gauss);
+                if(enthalpy>matpar->PureIceEnthalpy(pressure)){
+                        enthalpy_input->GetInputValue(&enthalpyup,gaussup);
+                        pressure_input->GetInputValue(&pressureup,gaussup);
+                        if(enthalpyup>matpar->PureIceEnthalpy(pressureup)){
+                                //do nothing, don't add heatflux
+                        }
+                        else{
+                                //need to change spcenthalpy according to Aschwanden 
+                                //NEED TO UPDATE
+                        }
+                }
+                else{
+                        geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
+                        friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
+                        vx_input->GetInputValue(&vx,gauss);
+                        vy_input->GetInputValue(&vy,gauss);
+                        basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
+
+                        scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
+                        if(dt) scalar=dt*scalar;
+
+                        for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
+                }
         }
 
         /*Clean up and return*/
         delete gauss;
+        delete gaussup;
         delete friction;
         return pe;
@@ -4085 +4107 @@
         double     rho_ice,heatcapacity,geothermalflux_value;
         double     basalfriction,alpha2,vx,vy;
-        double     scalar,enthalpy,enthalpyup;
-        double     pressure,pressureup;
         double     basis[NUMVERTICES];
+        double     scalar;
         Friction*  friction=NULL;
         GaussPenta* gauss=NULL;
-        GaussPenta* gaussup=NULL;
 
         /* Geothermal flux on ice sheet base and basal friction */
@@ -4108 +4128 @@
         Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
         Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
-        Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
-        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
         Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
 
@@ -4117 +4135 @@
         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
         gauss=new GaussPenta(0,1,2,2);
-        gaussup=new GaussPenta(3,4,5,2);
         for(ig=gauss->begin();ig<gauss->end();ig++){
 
                 gauss->GaussPoint(ig);
-                gaussup->GaussPoint(ig);
 
                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
                 GetNodalFunctionsP1(&basis[0], gauss);
 
-                enthalpy_input->GetInputValue(&enthalpy,gauss);
-                pressure_input->GetInputValue(&pressure,gauss);
-//              if(enthalpy>matpar->PureIceEnthalpy(pressure)){
-//                      enthalpy_input->GetInputValue(&enthalpyup,gaussup);
-//                      pressure_input->GetInputValue(&pressureup,gaussup);
-//                      if(enthalpyup>matpar->PureIceEnthalpy(pressureup)){
-//                              //do nothing, don't add heatflux
-//                      }
-//                      else{
-//                              //need to change spcenthalpy according to Aschwanden 
-//                              //NEED TO UPDATE
-//                      }
-//              }
-//              else{
                         geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
                         friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
@@ -4150 +4152 @@
 
                         for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
-//              }
         }
 
         /*Clean up and return*/
         delete gauss;
-        delete gaussup;
         delete friction;
         return pe;