Changeset 21406

Timestamp:

11/21/16 14:18:11 (8 years ago)

Author:

jbondzio

Message:

BUG: checking consistency of 'if waterfraction>0 then T=Tpmp', updating archives

Location:

issm/trunk-jpl

Files:

• src/m/classes/initialization.m (modified) (2 diffs)
• src/m/classes/initialization.py (modified) (2 diffs)
• test/Archives/Archive803.arch (modified) ( previous)
• test/Archives/Archive805.arch (modified) ( previous)

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

@@ -76 +76 @@
                                 md = checkfield(md,'fieldname','initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices 1]);
                                 md = checkfield(md,'fieldname','initialization.watercolumn'  ,'>=',0,'size',[md.mesh.numberofvertices 1]);
+                                pos=find(md.initialization.waterfraction>0.);
+                                if(~isempty(pos)),
+                                        md = checkfield(md,'fieldname', 'delta Tpmp', 'field',abs(md.initialization.temperature(pos)-(md.materials.meltingpoint-md.materials.beta*md.initialization.pressure(pos))),'<',1e-11,...
+                                        'message','set temperature to pressure melting point at locations with waterfraction>0');
+                                end
                         end
                         if ismember('HydrologyShreveAnalysis',analyses),
@@ -129 +134 @@
                         if md.thermal.isenthalpy,
                                 tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;
-                                pos  = find(md.initialization.temperature>=tpmp-1e-12);
+                                pos  = find(md.initialization.waterfraction>0.);
                                 enthalpy      = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature);
                                 enthalpy(pos) = md.materials.heatcapacity*(tpmp(pos) - md.constants.referencetemperature) + md.materials.latentheat*md.initialization.waterfraction(pos);

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

@@ -100 +100 @@
                                 md = checkfield(md,'fieldname','initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices])
                                 md = checkfield(md,'fieldname','initialization.watercolumn'  ,'>=',0,'size',[md.mesh.numberofvertices])
+                                pos = np.nonzero(md.initialization.waterfraction > 0.)[0]
+                                if(pos.size):
+                                        md = checkfield(md,'fieldname', 'delta Tpmp', 'field', np.absolute(md.initialization.temperature[pos]-(md.materials.meltingpoint-md.materials.beta*md.initialization.pressure[pos])),'<',1e-11, 'message','set temperature to pressure melting point at locations with waterfraction>0');
                 if 'HydrologyShreveAnalysis' in analyses:
                         if hasattr(md.hydrology,'hydrologyshreve'):
@@ -131 +134 @@
                 if md.thermal.isenthalpy:
                         tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;
-                        pos  = np.nonzero(md.initialization.temperature > tpmp)[0]
+                        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,)
+                        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');