Changeset 25171

Timestamp:

06/28/20 00:26:36 (5 years ago)

Author:

jdquinn

Message:

BUG: Fix for test2002; additional MATLAB -> Python translation; cleanup

Location:

issm/trunk-jpl

Files:

• src/m/classes/matenhancedice.py (modified) (7 diffs)
• src/m/classes/matestar.py (modified) (6 diffs)
• src/m/classes/rotational.py (modified) (1 diff)
• src/m/classes/slr.py (modified) (3 diffs)
• src/m/classes/surfaceload.py (modified) (1 diff)
• test/NightlyRun/test2002.py (modified) (1 diff)

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

@@ -9 +9 @@
     MATICE class definition
 
-       Usage:
-          matenhancedice = matenhancedice()
+        Usage:
+            matenhancedice = matenhancedice()
     """
 
-    def __init__(self):  # {{{
+    def __init__(self): #{{{
         self.rho_ice = 0.
         self.rho_water = 0.
@@ -32 +32 @@
         self.rheology_law = ''
 
-    #giaivins:
+        #GIA
         self.lithosphere_shear_modulus = 0.
         self.lithosphere_density = 0.
@@ -38 +38 @@
         self.mantle_density = 0.
 
-    #SLR
-        self.earth_density = 0  # average density of the Earth, (kg / m^3)
+        #SLR
+        self.earth_density = 0  # average density of the Earth, (kg/m^3)
         self.setdefaultparameters()
     #}}}
 
-    def __repr__(self):  # {{{
-        string = "   Materials:"
-        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_E", "enhancement factor"))
-        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' or 'LliboutryDuval'"))
-        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
+    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 [N s/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_E", "enhancement factor"))
+        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' or 'LliboutryDuval'"))
+        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]"))
+
+        return s
     #}}}
 
-    def extrude(self, md):  # {{{
+    def extrude(self, md): #{{{
         self.rheology_E = project3d(md, 'vector', self.rheology_E, 'type', 'node')
         self.rheology_B = project3d(md, 'vector', self.rheology_B, 'type', 'node')
@@ -77 +78 @@
     #}}}
 
-    def setdefaultparameters(self):  # {{{
+    def setdefaultparameters(self): #{{{
         #ice density (kg / m^3)
         self.rho_ice = 917.
@@ -120 +121 @@
     #}}}
 
-    def checkconsistency(self, md, solution, analyses):  # {{{
+    def checkconsistency(self, md, solution, analyses): #{{{
         md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
         md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0)
@@ -136 +137 @@
             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)
@@ -141 +143 @@
     # }}}
 
-    def marshall(self, prefix, md, fid):  # {{{
+    def marshall(self, prefix, md, fid): #{{{
         WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 4, 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double')

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

@@ +1 @@
+import numpy as np
+
+from checkfield import checkfield
 from fielddisplay import fielddisplay
 from project3d import project3d
-from checkfield import checkfield
 from WriteData import WriteData
 
 
 class matestar(object):
-    """
-    matestar class definition
+    '''
+    MATESTAR class definition
 
-       Usage:
-          matestar = matestar()
-    """
+        Usage:
+            matestar = matestar()
+    '''
 
-    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_Ec = float('NaN')
-        self.rheology_Es = float('NaN')
+        self.rheology_B = np.nan
+        self.rheology_Ec = np.nan
+        self.rheology_Es = np.nan
         self.rheology_law = ''
 
-        #giaivins:
+        #GIA
         self.lithosphere_shear_modulus = 0.
         self.lithosphere_density = 0.
@@ -38 +40 @@
         self.mantle_density = 0.
 
-        #slr
+        #SLR
         self.earth_density = 0
 
-        #set default parameters:
+        #set default parameters
         self.setdefaultparameters()
     #}}}
 
-    def __repr__(self):  # {{{
-        string = "   Materials:"
-        string = "%s\n%s" % (string, fielddisplay(self, 'rho_ice', 'ice density [kg / m^3]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'rho_water', 'ocean 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 / kg]'))
-        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 / 3)]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'rheology_Ec', 'compressive enhancement factor'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'rheology_Es', 'shear enhancement factor'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'rheology_law', ['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', ''Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']))
-        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
+    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', 'ocean 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 [N s/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 / kg]'))
+        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/3)]'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'rheology_Ec', 'compressive enhancement factor'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'rheology_Es', 'shear enhancement factor'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'rheology_law', ['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', ''Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']))
+        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]'))
+
+        return s
     #}}}
 
-    def extrude(self, md):  # {{{
+    def extrude(self, md): #{{{
         self.rheology_B = project3d(md, 'vector', self.rheology_B, 'type', 'node')
         self.rheology_Ec = project3d(md, 'vector', self.rheology_Ec, 'type', 'node')
         self.rheology_Es = project3d(md, 'vector', self.rheology_Es, 'type', 'node')
+
         return self
     #}}}
 
-    def setdefaultparameters(self):  # {{{
+    def setdefaultparameters(self): #{{{
         #ice density (kg / m^3)
         self.rho_ice = 917.
@@ -109 +113 @@
         #available: none, paterson and arrhenius
         self.rheology_law = 'Paterson'
-    # GIA:
+        # GIA
         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
+        #SLR
         self.earth_density = 5512  # average density of the Earth, (kg / m^3)
 
@@ -120 +124 @@
     #}}}
 
-    def checkconsistency(self, md, solution, analyses):  # {{{
+    def checkconsistency(self, md, solution, analyses): #{{{
         md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
         md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0)
@@ -143 +147 @@
     # }}}
 
-    def marshall(self, prefix, md, fid):  # {{{
+    def marshall(self, prefix, md, fid): #{{{
         WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 2, 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double')

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

@@ -46 +46 @@
 
     def checkconsistency(self, md, solution, analyses): # {{{
-        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isrotational == 0):
+        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0):
             return md
 

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

@@ -122 +122 @@
     def checkconsistency(self, md, solution, analyses):  # {{{
         #Early return
-        if (solution != 'SealevelriseAnalysis') or (solution == 'TransientSolution' and not md.transient.isslr):
+        if (solution != 'SealevelriseAnalysis') or (solution == 'TransientSolution' and md.transient.isslr == 0):
             return md
 
@@ -147 +147 @@
 
         #check that love numbers are provided at the same level of accuracy:
-        if (size(self.love_h, 0) != size(self.love_k, 0) | size(self.love_h, 0) != size(self.love_l, 0)):
-            error('slr error message: love numbers should be provided at the same level of accuracy')
+        if (size(self.love_h, 0) != size(self.love_k, 0)) or (size(self.love_h, 0) != size(self.love_l, 0)):
+            raise Exception('slr error message: love numbers should be provided at the same level of accuracy')
 
         #cross check that whereever we have an ice load, the mask is < 0 on each vertex:
@@ -160 +160 @@
         #a coupler to a planet model is provided.
         if self.geodetic and not md.transient.iscoupler and domaintype(md.mesh) != 'mesh3dsurface':
-            error('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!')
+            raise Exception('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!')
         return md
     # }}}

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

@@ -41 +41 @@
 
     def checkconsistency(self, md, solution, analyses): # {{{
-        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.issurfaceload == 0):
+        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0):
             return md
 
-        if len(self.icethicknesschange) > 0:
+        if len(self.icethicknesschange):
             md  = checkfield(md,'fieldname', 'solidearth.surfaceload.icethicknesschange', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
 
-        if len(self.waterheightchange) > 0:
+        if len(self.waterheightchange):
             md  = checkfield(md,'fieldname', 'solidearth.surfaceload.waterheightchange', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
 
-        if len(self.other) > 0:
+        if len(self.other):
             md  = checkfield(md,'fieldname', 'solidearth.surfaceload.other', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
 

issm/trunk-jpl/test/NightlyRun/test2002.py

@@ -20 +20 @@
 md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, 1))
 md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, 1))
-md.dsl.global_average_thermosteric_sea_level_change = np.zeros((1, 1))
+md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1))
 md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1))
 md.dsl.sea_water_pressure_change_at_sea_floor = np.zeros((md.mesh.numberofvertices + 1, 1))