Changeset 23812

Timestamp:

03/29/19 14:18:15 (6 years ago)

Author:

schlegel

Message:

CHG: add enthalpy to match the matlab version, for transient restart

File:

• issm/trunk-jpl/src/m/classes/initialization.py (modified) (4 diffs)

issm/trunk-jpl/src/m/classes/initialization.py

@@ -20 +20 @@
                 self.vz            = float('NaN')
                 self.vel           = float('NaN')
+                self.enthalpy      = float('NaN')
                 self.pressure      = float('NaN')
                 self.temperature   = float('NaN')
@@ -40 +41 @@
                 string="%s\n%s"%(string,fielddisplay(self,'pressure','pressure [Pa]'))
                 string="%s\n%s"%(string,fielddisplay(self,'temperature','temperature [K]'))
+                string="%s\n%s"%(string,fielddisplay(self,'enthalpy','enthalpy [J]'))
                 string="%s\n%s"%(string,fielddisplay(self,'waterfraction','fraction of water in the ice'))
                 string="%s\n%s"%(string,fielddisplay(self,'watercolumn','thickness of subglacial water [m]'))
@@ -54 +56 @@
                 self.vel=project3d(md,'vector',self.vel,'type','node')
                 self.temperature=project3d(md,'vector',self.temperature,'type','node')
+                self.enthalpy=project3d(md,'vector',self.enthalpy,'type','node')
                 self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node')
                 self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node')
@@ -133 +136 @@
                 
                 if md.thermal.isenthalpy:
+                    if (np.size(self.enthalpy)<=1):
                         tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;
                         pos  = np.nonzero(md.initialization.waterfraction > 0.)[0]
-                        enthalpy      = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature);
-                        enthalpy[pos] = md.materials.heatcapacity*(tpmp[pos].reshape(-1,) - md.constants.referencetemperature) + md.materials.latentheat*md.initialization.waterfraction[pos].reshape(-1,)
-                        WriteData(fid,prefix,'data',enthalpy,'format','DoubleMat','mattype',1,'name','md.initialization.enthalpy');
+                        self.enthalpy      = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature);
+                        self.enthalpy[pos] = md.materials.heatcapacity*(tpmp[pos].reshape(-1,) - md.constants.referencetemperature) + md.materials.latentheat*md.initialization.waterfraction[pos].reshape(-1,)
+
+                    WriteData(fid,prefix,'data',self.enthalpy,'format','DoubleMat','mattype',1,'name','md.initialization.enthalpy');
 
         # }}}