Changeset 20272

Timestamp:

02/29/16 10:47:12 (9 years ago)

Author:

jbondzio

Message:

BUG: CreatePVectorSheet uses enthalpy values interpolated at gauss-points, which may lead to inconsistency between applied basal boundary conditions. Keep applying basal heat flux at temperate base, since heated base becomes spc'd anyway.

File:

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

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

@@ -821 +821 @@
         if(!element->IsOnBase() || element->IsFloating()) return NULL;
 
-        bool isdynamicbasalspc;
+        bool converged, isdynamicbasalspc;
         int i, state;
+        int enthalpy_enum;
         IssmDouble  dt,Jdet,scalar;
         IssmDouble      enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
@@ -840 +841 @@
         element->FindParam(&dt,TimesteppingTimeStepEnum);
         element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
+        element->GetInputValue(&converged,ConvergedEnum);
+        if(dt==0. && !converged) enthalpy_enum=EnthalpyPicardEnum; // use enthalpy from last iteration
+        else enthalpy_enum=EnthalpyEnum; // use enthalpy from last time step
         Input* vx_input             = element->GetInput(VxEnum);                          _assert_(vx_input);
         Input* vy_input             = element->GetInput(VyEnum);                          _assert_(vy_input);
         Input* vz_input             = element->GetInput(VzEnum);                          _assert_(vz_input);
-        Input* enthalpy_input            = element->GetInput(EnthalpyPicardEnum);                                        _assert_(enthalpy_input);
+        Input* enthalpy_input            = element->GetInput(enthalpy_enum);                                     _assert_(enthalpy_input);
         Input* pressure_input            = element->GetInput(PressureEnum);                                                      _assert_(pressure_input);
         Input* watercolumn_input         = element->GetInput(WatercolumnEnum);                                                   _assert_(watercolumn_input);
@@ -876 +880 @@
 
                 switch (state) {
-                        case 0:
-                                // cold, dry base: apply basal surface forcing
+                        case 0: case 1: case 2: case 3:
+                                // cold, dry base; cold, wet base; refreezing temperate base; thin temperate base: 
+                                // Apply basal surface forcing.
+                                // Interpolated values of enthalpy on gauss nodes may indicate cold base, 
+                                // although one node might have become temperate. So keep heat flux switched on.
                                 geothermalflux_input->GetInputValue(&geothermalflux,gauss);
                                 friction->GetAlpha2(&alpha2,gauss);
@@ -890 +897 @@
                                         pe->values[i]+=scalar*basis[i];
                                 break;
-                        case 1:
-                                // cold, wet base: keep at pressure melting point 
-                                break;
-                        case 2:
-                                // temperate, thin refreezing base: release spc
-                                break;
-                        case 3:
-                                // temperate, thin melting base: set spc
-                                break;
                         case 4:
                                 // temperate, thick melting base: set grad H*n=0