Changeset 25836 for issm/trunk/src/m/classes/matice.py

Timestamp:

12/08/20 08:45:53 (4 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 25834

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/m/classes/matice.py (modified) (8 diffs)

issm/trunk/src/m/classes/matice.py

@@ +1 @@
+import numpy as np
+
 from fielddisplay import fielddisplay
 from project3d import project3d
@@ -6 +8 @@
 
 class matice(object):
-    """
+    '''
     MATICE class definition
 
-       Usage:
-          matice = matice()
-    """
+    Usage:
+            matice = matice()
+    '''
 
-    def __init__(self):  # {{{
+    def __init__(self): #{{{
         self.rho_ice = 0.
         self.rho_water = 0.
@@ -27 +29 @@
         self.mixed_layer_capacity = 0.
         self.thermal_exchange_velocity = 0.
-        self.rheology_B = float('NaN')
-        self.rheology_n = float('NaN')
+        self.rheology_B = np.nan
+        self.rheology_n = np.nan
         self.rheology_law = ''
 
-        #giaivins:
+        #giaivins
         self.lithosphere_shear_modulus = 0.
         self.lithosphere_density = 0.
@@ -37 +39 @@
         self.mantle_density = 0.
 
-        #SLR
-        self.earth_density = 5512
+        #slr
+        self.earth_density = 0
         self.setdefaultparameters()
     #}}}
 
-    def __repr__(self):  # {{{
-        string = "   Materials:"
+    def __repr__(self): #{{{
+        s = "   Materials:"
+        s = "%s\n%s" % (s, fielddisplay(self, "rho_ice", "ice density [kg/m^3]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "rho_water", "water density [kg/m^3]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "rho_freshwater", "fresh water density [kg/m^3]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "mu_water", "water viscosity [Ns/m^2]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "heatcapacity", "heat capacity [J/kg/K]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "thermalconductivity", "ice thermal conductivity [W/m/K]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "temperateiceconductivity", "temperate ice thermal conductivity [W/m/K]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "effectiveconductivity_averaging", "computation of effectiveconductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)"))
+        s = "%s\n%s" % (s, fielddisplay(self, "meltingpoint", "melting point of ice at 1atm in K"))
+        s = "%s\n%s" % (s, fielddisplay(self, "latentheat", "latent heat of fusion [J/m^3]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "beta", "rate of change of melting point with pressure [K/Pa]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "mixed_layer_capacity", "mixed layer capacity [W/kg/K]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "thermal_exchange_velocity", "thermal exchange velocity [m/s]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1/n)]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
+        s = "%s\n%s" % (s, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', or 'NyeH2O'"))
+        s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_density", "Lithosphere density [g/cm^-3]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "mantle_density", "Mantle density [g/cm^-3]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "earth_density", "Mantle density [kg/m^-3]"))
 
-        string = "%s\n%s" % (string, fielddisplay(self, "rho_ice", "ice density [kg / m^3]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "rho_water", "water density [kg / m^3]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "rho_freshwater", "fresh water density [kg / m^3]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "mu_water", "water viscosity [N s / m^2]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "heatcapacity", "heat capacity [J / kg / K]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "thermalconductivity", "ice thermal conductivity [W / m / K]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "temperateiceconductivity", "temperate ice thermal conductivity [W / m / K]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "effectiveconductivity_averaging", "computation of effectiveconductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)"))
-        string = "%s\n%s" % (string, fielddisplay(self, "meltingpoint", "melting point of ice at 1atm in K"))
-        string = "%s\n%s" % (string, fielddisplay(self, "latentheat", "latent heat of fusion [J / m^3]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "beta", "rate of change of melting point with pressure [K / Pa]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "mixed_layer_capacity", "mixed layer capacity [W / kg / K]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "thermal_exchange_velocity", "thermal exchange velocity [m / s]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1 / n)]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
-        string = "%s\n%s" % (string, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', or 'NyeH2O'"))
-        string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_density", "Lithosphere density [g / cm^ - 3]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "mantle_density", "Mantle density [g / cm^ - 3]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "earth_density", "Mantle density [kg / m^ - 3]"))
-        return string
+        return s
     #}}}
 
-    def extrude(self, md):  # {{{
+    def extrude(self, md): #{{{
         self.rheology_B = project3d(md, 'vector', self.rheology_B, 'type', 'node')
         self.rheology_n = project3d(md, 'vector', self.rheology_n, 'type', 'element')
@@ -75 +77 @@
     #}}}
 
-    def setdefaultparameters(self):  # {{{
-        #ice density (kg / m^3)
+    def setdefaultparameters(self): #{{{
+        #ice density (kg/m^3)
         self.rho_ice = 917.
-        #ocean water density (kg / m^3)
+        #ocean water density (kg/m^3)
         self.rho_water = 1023.
-        #fresh water density (kg / m^3)
+        #fresh water density (kg/m^3)
         self.rho_freshwater = 1000.
-        #water viscosity (N.s / m^2)
+        #water viscosity (N.s/m^2)
         self.mu_water = 0.001787
-        #ice heat capacity cp (J / kg / K)
+        #ice heat capacity cp (J/kg/K)
         self.heatcapacity = 2093.
-        #ice latent heat of fusion L (J / kg)
-        self.latentheat = 3.34 * 1.0e5
-        #ice thermal conductivity (W / m / K)
+        #ice latent heat of fusion L (J/kg)
+        self.latentheat = 3.34e5
+        #ice thermal conductivity (W/m/K)
         self.thermalconductivity = 2.4
+        #temperate ice thermal conductivity (W/m/K)
+        self.temperateiceconductivity = 0.24
         #computation of effective conductivity
         self.effectiveconductivity_averaging = 1
-        #temperate ice thermal conductivity (W / m / K)
-        self.temperateiceconductivity = 0.24
         #the melting point of ice at 1 atmosphere of pressure in K
         self.meltingpoint = 273.15
-        #rate of change of melting point with pressure (K / Pa)
-        self.beta = 9.8 * 1.0e-8
-        #mixed layer (ice-water interface) heat capacity (J / kg / K)
+        #rate of change of melting point with pressure (K/Pa)
+        self.beta = 9.8e-8
+        #mixed layer (ice-water interface) heat capacity (J/kg/K)
         self.mixed_layer_capacity = 3974.
-        #thermal exchange velocity (ice-water interface) (m / s)
-        self.thermal_exchange_velocity = 1.00 * 1.0e-4
+        #thermal exchange velocity (ice-water interface) (m/s)
+        self.thermal_exchange_velocity = 1.00e-4
         #Rheology law: what is the temperature dependence of B with T
         #available: none, paterson and arrhenius
@@ -107 +109 @@
 
         # GIA:
-        self.lithosphere_shear_modulus = 6.7 * 1.0e10  # (Pa)
-        self.lithosphere_density = 3.32  # (g / cm^ - 3)
-        self.mantle_shear_modulus = 1.45 * 1.0e11  # (Pa)
-        self.mantle_density = 3.34  # (g / cm^ - 3)
+        self.lithosphere_shear_modulus = 6.7e10 # (Pa)
+        self.lithosphere_density = 3.32  # (g/cm^-3)
+        self.mantle_shear_modulus = 1.45e11  # (Pa)
+        self.mantle_density = 3.34  # (g/cm^-3)
 
-        #SLR
-        self.earth_density = 5512  # average density of the Earth, (kg / m^3)
-        return self
+        # SLR
+        self.earth_density = 5512  # average density of the Earth, (kg/m^3)
     #}}}
 
-    def checkconsistency(self, md, solution, analyses):  # {{{
-        md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
-        md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0)
-        md = checkfield(md, 'fieldname', 'materials.rho_freshwater', '>', 0)
-        md = checkfield(md, 'fieldname', 'materials.mu_water', '>', 0)
-        md = checkfield(md, 'fieldname', 'materials.rheology_B', '>', 0, 'timeseries', 1, 'NaN', 1, 'Inf', 1)
-        md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements])
-        md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O'])
-        md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
-        md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', [1])
-        md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', [1])
-        md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', [1])
-        md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
-        md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
+    def checkconsistency(self, md, solution, analyses): #{{{
+        if solution != 'SealevelriseSolution':
+            md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
+            md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0)
+            md = checkfield(md, 'fieldname', 'materials.rho_freshwater', '>', 0)
+            md = checkfield(md, 'fieldname', 'materials.mu_water', '>', 0)
+            md = checkfield(md, 'fieldname', 'materials.rheology_B', '>', 0, 'universal', 1, 'NaN', 1, 'Inf', 1)
+            md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements])
+            md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O'])
+            md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
+
+        if 'GiaAnalysis' in analyses:
+            md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', [1])
+            md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', [1])
+            md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', [1])
+            md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
+
+        if 'SealevelriseAnalysis' in analyses:
+            md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
+
         return md
-    # }}}
+    #}}}
 
-    def marshall(self, prefix, md, fid):  # {{{
+    def marshall(self, prefix, md, fid): #{{{
         WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 3, 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double')
@@ -152 +159 @@
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2)
         WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
-
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10.**3.)
@@ -158 +164 @@
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10.**3.)
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
-
-    # }}}
+    #}}}