Changeset 16621 for issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

Timestamp:

11/06/13 10:02:04 (12 years ago)

Author:

jbondzio

Message:

CHG: Added harmonic mean for better convergence of steadystate and transient solver when handling a jump in enthalpy conductivity

File:

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

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

@@ -3776 +3776 @@
         IssmDouble epsvel=2.220446049250313e-16;
         IssmDouble heatcapacity,thermalconductivity,dt;
-        IssmDouble pressure[NUMVERTICES],enthalpy[NUMVERTICES];
+        IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES];
         IssmDouble latentheat,kappa;
         IssmDouble tau_parameter,diameter;
@@ -3824 +3824 @@
                 /*Conduction: */  
                 GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); 
-
-                //enthalpy_input->GetInputValue(&enthalpy, gauss);
-                //pressure_input->GetInputValue(&pressure, gauss);
-                GetInputListOnVertices(&enthalpy[0],EnthalpyPicardEnum);
+                GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
                 GetInputListOnVertices(&pressure[0],PressureEnum);
-                kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy,pressure); _assert_(kappa>0.);  // TODO: switch to pointers
+                kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy,pressure); _assert_(kappa>0.);
 
                 D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice;
@@ -4232 +4229 @@
         IssmDouble rho_ice,heatcapacity;
         IssmDouble thermalconductivity,kappa;
-        IssmDouble viscosity,pressure;
-        IssmDouble enthalpy;
-        IssmDouble enthalpy_v[NUMVERTICES], pressure_v[NUMVERTICES];
+        IssmDouble viscosity,pressure,enthalpy;
+        IssmDouble enthalpypicard[NUMVERTICES], pressure_p[NUMVERTICES];
         IssmDouble tau_parameter,diameter;
         IssmDouble u,v,w;
@@ -4300 +4296 @@
                         vy_input->GetInputValue(&v, gauss);
                         vz_input->GetInputValue(&w, gauss);
-                        //pressure_input->GetInputValue(&pressure, gauss);
-                        //enthalpypicard_input->GetInputValue(&enthalpypicard, gauss);
-                        GetInputListOnVertices(&enthalpy_v[0],EnthalpyPicardEnum);
-                        GetInputListOnVertices(&pressure_v[0],PressureEnum);
-                        kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy_v,pressure_v);
+                        GetInputListOnVertices(&enthalpypicard[0],EnthalpyPicardEnum);
+                        GetInputListOnVertices(&pressure_p[0],PressureEnum);
+                        kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpypicard,pressure_p);
                         kappa/=rho_ice;
                         tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);