Ignore:
Timestamp:
12/22/21 10:39:44 (3 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 26742

Location:
issm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/m/classes/materials.py

    r25836 r26744  
    88
    99class materials(object):
    10     '''
    11     MATERIALS class definition
    12 
    13        Usage:
    14           materials = materials()
    15     '''
     10    """MATERIALS class definition
     11
     12    Usage:
     13        materials = materials()
     14    """
    1615
    1716    def __init__(self, *args): #{{{
    1817        self.nature = []
    19         if not len(args):
     18        if len(args) == 0:
    2019            self.nature = ['ice']
    2120        else:
    … …  
    2625                raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
    2726
    28         #start filling in the dynamic fields:
    29         for i in range(len(self.nature)):
    30             nat = self.nature[i]
    31             if nat == 'ice':
    32                 setattr(self, 'rho_ice', 0)
    33                 setattr(self, 'rho_water', 0)
    34                 setattr(self, 'rho_freshwater', 0)
    35                 setattr(self, 'mu_water', 0)
    36                 setattr(self, 'heatcapacity', 0)
    37                 setattr(self, 'latentheat', 0)
    38                 setattr(self, 'thermalconductivity', 0)
    39                 setattr(self, 'temperateiceconductivity', 0)
    40                 setattr(self, 'meltingpoint', 0)
    41                 setattr(self, 'beta', 0)
    42                 setattr(self, 'mixed_layer_capacity', 0)
    43                 setattr(self, 'thermal_exchange_velocity', 0)
    44                 setattr(self, 'rheology_B', 0)
    45                 setattr(self, 'rheology_n', 0)
    46                 setattr(self, 'rheology_law', 0)
    47             elif nat == 'litho':
    48                 setattr(self, 'numlayers', 0)
    49                 setattr(self, 'radius', 0)
    50                 setattr(self, 'viscosity', 0)
    51                 setattr(self, 'lame_lambda', 0)
    52                 setattr(self, 'lame_mu', 0)
    53                 setattr(self, 'burgers_viscosity', 0)
    54                 setattr(self, 'burgers_mu', 0)
    55                 setattr(self, 'isburgers', 0)
    56                 setattr(self, 'density', 0)
    57                 setattr(self, 'issolid', 0)
    58             elif nat == 'hydro':
    59                 setattr(self, 'rho_ice', 0)
    60                 setattr(self, 'rho_water', 0)
    61                 setattr(self, 'rho_freshwater', 0)
    62                 setattr(self, 'earth_density', 0)
     27        # Start filling in the dynamic fields (not truly dynamic under Python)
     28        for i in range(len(self.nature)):
     29            nat = self.nature[i]
     30            if nat == 'ice':
     31                self.rho_ice = 0
     32                self.rho_water = 0
     33                self.rho_freshwater = 0
     34                self.mu_water = 0
     35                self.heatcapacity = 0
     36                self.latentheat = 0
     37                self.thermalconductivity = 0
     38                self.temperateiceconductivity = 0
     39                self.effectiveconductivity_averaging = 0
     40                self.meltingpoint = 0
     41                self.beta = 0
     42                self.mixed_layer_capacity = 0
     43                self.thermal_exchange_velocity = 0
     44                self.rheology_B = 0
     45                self.rheology_n = 0
     46                self.rheology_law = 0
     47            elif nat == 'litho':
     48                self.numlayers = 0
     49                self.radius = 0
     50                self.viscosity = 0
     51                self.lame_lambda = 0
     52                self.lame_mu = 0
     53                self.burgers_viscosity = 0
     54                self.burgers_mu = 0
     55                self.ebm_alpha = 0
     56                self.ebm_delta = 0
     57                self.ebm_taul = 0
     58                self.ebm_tauh = 0
     59                self.rheologymodel = 0
     60                self.density = 0
     61                self.issolid = 0
     62            elif nat == 'hydro':
     63                self.rho_ice = 0
     64                self.rho_water = 0
     65                self.rho_freshwater = 0
    6366            else:
    6467                raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
    65 
    66         #set default parameters:
     68        self.earth_density = 0
     69
    6770        self.setdefaultparameters()
    6871    #}}}
    6972
     73    def __repr__(self): #{{{
     74        s = '   Materials:\n'
     75        for i in range(len(self.nature)):
     76            nat = self.nature[i]
     77            if nat == 'ice':
     78                s += '\n      Ice:\n'
     79                s += '{}\n'.format(fielddisplay(self, 'rho_ice', 'ice density [kg/m^3]'))
     80                s += '{}\n'.format(fielddisplay(self, 'rho_water', 'ocean water density [kg/m^3]'))
     81                s += '{}\n'.format(fielddisplay(self, 'rho_freshwater', 'fresh water density [kg/m^3]'))
     82                s += '{}\n'.format(fielddisplay(self, 'mu_water', 'water viscosity [N s/m^2]'))
     83                s += '{}\n'.format(fielddisplay(self, 'heatcapacity', 'heat capacity [J/kg/K]'))
     84                s += '{}\n'.format(fielddisplay(self, 'thermalconductivity', 'ice thermal conductivity [W/m/K]'))
     85                s += '{}\n'.format(fielddisplay(self, 'temperateiceconductivity', 'temperate ice thermal conductivity [W/m/K]'))
     86                s += '{}\n'.format(fielddisplay(self, 'meltingpoint', 'melting point of ice at 1atm in K'))
     87                s += '{}\n'.format(fielddisplay(self, 'latentheat', 'latent heat of fusion [J/m^3]'))
     88                s += '{}\n'.format(fielddisplay(self, 'beta', 'rate of change of melting point with pressure [K/Pa]'))
     89                s += '{}\n'.format(fielddisplay(self, 'mixed_layer_capacity', 'mixed layer capacity [W/kg/K]'))
     90                s += '{}\n'.format(fielddisplay(self, 'thermal_exchange_velocity', 'thermal exchange velocity [m/s]'))
     91                s += '{}\n'.format(fielddisplay(self, 'rheology_B', 'flow law parameter [Pa s^(1/n)]'))
     92                s += '{}\n'.format(fielddisplay(self, 'rheology_n', 'Glen\'s flow law exponent'))
     93                s += '{}\n'.format(fielddisplay(self, 'rheology_law', 'law for the temperature dependance of the rheology: \'None\', \'BuddJacka\', \'Cuffey\', \'CuffeyTemperate\', \'Paterson\', \'Arrhenius\', \'LliboutryDuval\', \'NyeCO2\', or \'NyeH2O\''))
     94            elif nat == 'litho':
     95                s += '\n      Litho:\n'
     96                s += '{}\n'.format(fielddisplay(self, 'numlayers', 'number of layers (default: 2)'))
     97                s += '{}\n'.format(fielddisplay(self, 'radius', 'array describing the radius for each interface (numlayers + 1) [m]'))
     98                s += '{}\n'.format(fielddisplay(self, 'viscosity', 'array describing each layer\'s viscosity (numlayers) [Pa.s]'))
     99                s += '{}\n'.format(fielddisplay(self, 'lame_lambda', 'array describing the lame lambda parameter (numlayers) [Pa]'))
     100                s += '{}\n'.format(fielddisplay(self, 'lame_mu', 'array describing the shear modulus for each layers (numlayers) [Pa]'))
     101                s += '{}\n'.format(fielddisplay(self, 'burgers_viscosity', 'array describing each layer\'s transient viscosity, only for Burgers rheologies  (numlayers) [Pa.s]'))
     102                s += '{}\n'.format(fielddisplay(self, 'burgers_mu', 'array describing each layer\'s transient shear modulus, only for Burgers rheologies  (numlayers) [Pa]'))
     103
     104                s += '{}\n'.format(fielddisplay(self, 'ebm_alpha', 'array describing each layer\'s exponent parameter controlling the shape of shear modulus curve between taul and tauh, only for EBM rheology (numlayers)'))
     105                s += '{}\n'.format(fielddisplay(self, 'ebm_delta', 'array describing each layer\'s amplitude of the transient relaxation (ratio between elastic rigity to pre-maxwell relaxation rigity), only for EBM rheology (numlayers)'))
     106                s += '{}\n'.format(fielddisplay(self, 'ebm_taul', 'array describing each layer\'s starting period for transient relaxation, only for EBM rheology  (numlayers) [s]'))
     107                s += '{}\n'.format(fielddisplay(self, 'ebm_tauh', 'array describing each layer''s array describing each layer\'s end period for transient relaxation, only for Burgers rheology (numlayers) [s]'))
     108
     109                s += '{}\n'.format(fielddisplay(self, 'rheologymodel', 'array describing whether we adopt a Maxwell (0), Burgers (1) or EBM (2) rheology (default: 0)'))
     110                s += '{}\n'.format(fielddisplay(self, 'density', 'array describing each layer\'s density (numlayers) [kg/m^3]'))
     111                s += '{}\n'.format(fielddisplay(self, 'issolid', 'array describing whether the layer is solid or liquid (default: 1) (numlayers)'))
     112            elif nat == 'hydro':
     113                s += '\n      Hydro:\n'
     114                s += '{}\n'.format(fielddisplay(self, 'rho_ice', 'ice density [kg/m^3]'))
     115                s += '{}\n'.format(fielddisplay(self, 'rho_water', 'ocean water density [kg/m^3]'))
     116                s += '{}\n'.format(fielddisplay(self, 'earth_density', 'mantle density [kg/m^3]'))
     117                s += '{}\n'.format(fielddisplay(self, 'rho_freshwater', 'fresh water density [kg/m^3]'))
     118
     119            else:
     120                raise RuntimeError('materials constructor error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')')
     121        return s
     122    #}}}
     123
    70124    def setdefaultparameters(self): #{{{
    71125        for i in range(len(self.nature)):
    72126            nat = self.nature[i]
    73127            if nat == 'ice':
    74                 #ice density (kg/m^3)
    75                 self.rho_ice = 917.
    76                 #ocean water density (kg/m^3)
    77                 self.rho_water = 1023.
    78                 #fresh water density (kg/m^3)
    79                 self.rho_freshwater = 1000.
    80                 #water viscosity (N.s/m^2)
     128                # Ice density (kg/m^3)
     129                self.rho_ice = 917.0
     130
     131                # Ocean water density (kg/m^3)
     132                self.rho_water = 1023.0
     133
     134                # Fresh water density (kg/m^3)
     135                self.rho_freshwater = 1000.0
     136
     137                # Water viscosity (N.s/m^2)
    81138                self.mu_water = 0.001787
    82                 #ice heat capacity cp (J/kg/K)
    83                 self.heatcapacity = 2093.
    84                 #ice latent heat of fusion L (J / kg)
    85                 self.latentheat = 3.34e5
    86                 #ice thermal conductivity (W/m/K)
     139
     140                # Ice heat capacity cp (J/kg/K)
     141                self.heatcapacity = 2093.0
     142
     143                # Ice latent heat of fusion L (J/kg)
     144                self.latentheat = 3.34 * 1e5
     145
     146                # Ice thermal conductivity (W/m/K)
    87147                self.thermalconductivity = 2.4
    88                 #wet ice thermal conductivity (W/m/K)
     148
     149                # Wet ice thermal conductivity (W/m/K)
    89150                self.temperateiceconductivity = 0.24
    90                 #the melting point of ice at 1 atmosphere of pressure in K
     151
     152                # Computation of effective conductivity
     153                self.effectiveconductivity_averaging = 1
     154
     155                # The melting point of ice at 1 atmosphere of pressure in K
    91156                self.meltingpoint = 273.15
    92                 #rate of change of melting point with pressure (K/Pa)
    93                 self.beta = 9.8e-8
    94                 #mixed layer (ice-water interface) heat capacity (J/kg/K)
    95                 self.mixed_layer_capacity = 3974.
    96                 #thermal exchange velocity (ice-water interface) (m/s)
    97                 self.thermal_exchange_velocity = 1.00e-4
    98                 #Rheology law: what is the temperature dependence of B with T
    99                 #available: none, paterson and arrhenius
     157
     158                # Rate of change of melting point with pressure (K/Pa)
     159                self.beta = 9.8 * 1e-8
     160
     161                # Mixed layer (ice-water interface) heat capacity (J/kg/K)
     162                self.mixed_layer_capacity = 3974.0
     163
     164                # Thermal exchange velocity (ice-water interface) (m/s)
     165                self.thermal_exchange_velocity = 1.00 * 1e-4
     166
     167                # Rheology law: what is the temperature dependence of B with T
     168                # available: none, paterson and arrhenius
    100169                self.rheology_law = 'Paterson'
    101             elif nat == 'litho':
    102                 #we default to a configuration that enables running GIA solutions using giacaron and / or giaivins.
     170
     171                # Rheology fields default
     172                self.rheology_B = 1 * 1e8
     173                self.rheology_n = 3
     174            elif nat == 'litho':
     175                # We default to a configuration that enables running GIA
     176                # solutions using giacaron and/or giaivins
    103177                self.numlayers = 2
    104                 #center of the earth (approximation, must not be 0), then the lab (lithosphere / asthenosphere boundary) then the surface
    105                 #(with 1d3 to avoid numerical singularities)
    106                 self.radius = [1e3, 6278e3, 6378e3]
    107                 self.viscosity = [1e21, 1e40]  #mantle and lithosphere viscosity (respectively) [Pa.s]
    108                 self.lame_mu = [1.45e11, 6.7e10]  # (Pa)  #lithosphere and mantle shear modulus (respectively) [Pa]
    109                 self.lame_lambda = self.lame_mu  # (Pa)  #mantle and lithosphere lamba parameter (respectively) [Pa]
     178
     179                # Center of the earth (approximation, must not be 0), then the
     180                # lab (lithosphere/asthenosphere boundary) then the surface
     181                # (with 1d3 to avoid numerical singularities)
     182                self.radius = [1e3, 6278 * 1e3, 6378 * 1e3]
     183
     184                self.viscosity = [1e21, 1e40] # Mantle and lithosphere viscosity (respectively) [Pa.s]
     185                self.lame_mu = [1.45 * 1e11, 6.7 * 1e10] # (Pa) # Lithosphere and mantle shear modulus (respectively) [Pa]
     186                self.lame_lambda = self.lame_mu # (Pa) # Mantle and lithosphere lamba parameter (respectively) [Pa]
    110187                self.burgers_viscosity = [np.nan, np.nan]
    111188                self.burgers_mu = [np.nan, np.nan]
    112                 self.isburgers = [0, 0]
    113                 self.density = [5.51e3, 5.50e3]  # (Pa)  #mantle and lithosphere density [kg/m^3]
    114                 self.issolid = [1, 1]  # is layer solid or liquid.
    115             elif nat == 'hydro':
    116                 #ice density (kg/m^3)
    117                 self.rho_ice = 917.
    118                 #ocean water density (kg/m^3)
    119                 self.rho_water = 1023.
    120                 #average density of the Earth (kg/m^3)
    121                 self.earth_density = 5512
    122                 #fresh water density (kg/m^3)
    123                 self.rho_freshwater = 1000.
     189
     190                self.ebm_alpha = [np.nan, np.nan]
     191                self.ebm_delta = [np.nan, np.nan]
     192                self.ebm_taul = [np.nan, np.nan]
     193                self.ebm_tauh = [np.nan, np.nan]
     194                self.rheologymodel = [0, 0]
     195                self.density = [5.51 * 1e3, 5.50 * 1e3] # (Pa) # Mantle and lithosphere density [kg/m^3]
     196                self.issolid = [1, 1] # Is layer solid or liquid?
     197            elif nat == 'hydro':
     198                # Ice density (kg/m^3)
     199                self.rho_ice = 917.0
     200
     201                # Ocean water density (kg/m^3)
     202                self.rho_water = 1023.0
     203
     204                # Fresh water density (kg/m^3)
     205                self.rho_freshwater = 1000.0
    124206            else:
    125207                raise RuntimeError("materials setdefaultparameters error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
    126     #}}}
    127 
    128     def __repr__(self): #{{{
    129         string = "   Materials:"
    130         for i in range(len(self.nature)):
    131             nat = self.nature[i]
    132             if nat == 'ice':
    133                 string = "%s\n%s" % (string, 'Ice:')
    134                 string = "%s\n%s" % (string, fielddisplay(self, "rho_ice", "ice density [kg/m^3]"))
    135                 string = "%s\n%s" % (string, fielddisplay(self, "rho_water", "ocean water density [kg/m^3]"))
    136                 string = "%s\n%s" % (string, fielddisplay(self, "rho_freshwater", "fresh water density [kg/m^3]"))
    137                 string = "%s\n%s" % (string, fielddisplay(self, "mu_water", "water viscosity [N s/m^2]"))
    138                 string = "%s\n%s" % (string, fielddisplay(self, "heatcapacity", "heat capacity [J/kg/K]"))
    139                 string = "%s\n%s" % (string, fielddisplay(self, "thermalconductivity", "ice thermal conductivity [W/m/K]"))
    140                 string = "%s\n%s" % (string, fielddisplay(self, "temperateiceconductivity", "temperate ice thermal conductivity [W/m/K]"))
    141                 string = "%s\n%s" % (string, fielddisplay(self, "meltingpoint", "melting point of ice at 1atm in K"))
    142                 string = "%s\n%s" % (string, fielddisplay(self, "latentheat", "latent heat of fusion [J/m^3]"))
    143                 string = "%s\n%s" % (string, fielddisplay(self, "beta", "rate of change of melting point with pressure [K/Pa]"))
    144                 string = "%s\n%s" % (string, fielddisplay(self, "mixed_layer_capacity", "mixed layer capacity [W/kg/K]"))
    145                 string = "%s\n%s" % (string, fielddisplay(self, "thermal_exchange_velocity", "thermal exchange velocity [m/s]"))
    146                 string = "%s\n%s" % (string, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1/n)]"))
    147                 string = "%s\n%s" % (string, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
    148                 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'"))
    149             elif nat == 'litho':
    150                 string = "%s\n%s" % (string, 'Litho:')
    151                 string = "%s\n%s" % (string, fielddisplay(self, 'numlayers', 'number of layers (default 2)'))
    152                 string = "%s\n%s" % (string, fielddisplay(self, 'radius', 'array describing the radius for each interface (numlayers + 1) [m]'))
    153                 string = "%s\n%s" % (string, fielddisplay(self, 'viscosity', 'array describing each layer''s viscosity (numlayers) [Pa.s]'))
    154                 string = "%s\n%s" % (string, fielddisplay(self, 'lame_lambda', 'array describing the lame lambda parameter (numlayers) [Pa]'))
    155                 string = "%s\n%s" % (string, fielddisplay(self, 'lame_mu', 'array describing the shear modulus for each layers (numlayers) [Pa]'))
    156                 string = "%s\n%s" % (string, fielddisplay(self, 'burgers_viscosity', 'array describing each layer''s transient viscosity, only for Burgers rheologies  (numlayers) [Pa.s]'))
    157                 string = "%s\n%s" % (string, fielddisplay(self, 'burgers_mu', 'array describing each layer''s transient shear modulus, only for Burgers rheologies  (numlayers) [Pa]'))
    158                 string = "%s\n%s" % (string, fielddisplay(self, 'isburgers', 'array describing whether we adopt a MaxWell (0) or Burgers (1) rheology (default 0)'))
    159                 string = "%s\n%s" % (string, fielddisplay(self, 'density', 'array describing each layer''s density (numlayers) [kg/m^3]'))
    160                 string = "%s\n%s" % (string, fielddisplay(self, 'issolid', 'array describing whether the layer is solid or liquid (default 1) (numlayers)'))
    161             elif nat == 'hydro':
    162                 string = "%s\n%s" % (string, 'Hydro:')
    163                 string = "%s\n%s" % (string, fielddisplay(self, "rho_ice", "ice density [kg/m^3]"))
    164                 string = "%s\n%s" % (string, fielddisplay(self, "rho_water", "ocean water density [kg/m^3]"))
    165                 string = "%s\n%s" % (string, fielddisplay(self, "earth_density", "mantle density [kg/m^3]"))
    166                 string = "%s\n%s" % (string, fielddisplay(self, "rho_freshwater", "fresh water density [kg/m^3]"))
    167 
    168             else:
    169                 raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
    170 
    171         return string
     208
     209            # Average density of the Earth (kg/m^3)
     210            self.earth_density = 5512
    172211    #}}}
    173212
    … …  
    194233                md = checkfield(md, 'fieldname', 'materials.density', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>', 0)
    195234                md = checkfield(md, 'fieldname', 'materials.viscosity', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
    196                 md = checkfield(md, 'fieldname', 'materials.isburgers', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0, '<=', 1)
     235                md = checkfield(md, 'fieldname', 'materials.rheologymodel', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0, '<=', 2)
    197236                md = checkfield(md, 'fieldname', 'materials.burgers_viscosity', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
    198237                md = checkfield(md, 'fieldname', 'materials.burgers_mu', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
     238                md = checkfield(md, 'fieldname', 'materials.ebm_alpha', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
     239                md = checkfield(md, 'fieldname', 'materials.ebm_delta', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
     240                md = checkfield(md, 'fieldname', 'materials.ebm_taul', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
     241                md = checkfield(md, 'fieldname', 'materials.ebm_tauh', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
    199242
    200243                for i in range(md.materials.numlayers):
    201                     if md.materials.isburgers[i] and (np.isnan(md.materials.burgers_viscosity[i] or np.isnan(md.materials.burgers_mu[i]))):
    202                         raise RuntimeError("materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with isburgers choice")
    203 
    204                     if md.materials.issolid[0] == 0 or md.materials.lame_mu[0] == 0:
    205                         raise RuntimeError('First layer must be solid (issolid(1) > 0 AND lame_mu(1) > 0). Add a weak inner core if necessary.')
    206 
    207                     for i in range(md.materials.numlayers - 1):
    208                         if (not md.materials.issolid[i]) and (not md.materials.issolid[i + 1]):  #if there are at least two consecutive indices that contain issolid = 0
    209                             raise RuntimeError("%s%i%s" % ('2 or more adjacent fluid layers detected starting at layer ', i, '. This is not supported yet. Consider merging them.'))
     244                    if md.materials.rheologymodel[i] == 1 and (np.isnan(md.materials.burgers_viscosity[i] or np.isnan(md.materials.burgers_mu[i]))):
     245                        raise RuntimeError('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with rheologymodel choice')
     246
     247                    if md.materials.rheologymodel[i] == 2 and (np.isnan(md.materials.ebm_alpha[i]) or np.isnan(md.materials.ebm_delta[i]) or np.isnan(md.materials.ebm_taul[i]) or np.isnan(md.materials.ebm_tauh[i])):
     248                        raise RuntimeError('materials checkconsistency error message: Litho ebm_alpha, ebm_delta, ebm_taul or ebm_tauh has NaN values, inconsistent with rheologymodel choice')
     249                if md.materials.issolid[0] == 0 or md.materials.lame_mu[0] == 0:
     250                    raise RuntimeError('First layer must be solid (issolid[0] > 0 AND lame_mu[0] > 0). Add a weak inner core if necessary.')
     251                ind = np.where(md.materials.issolid == 0)[0]
     252                if np.sum(np.in1d(np.diff(ind),1) >= 1): # If there are at least two consecutive indices that contain issolid = 0
     253                    raise RuntimeError('Fluid layers detected at layers #' + indices + ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.')
     254
    210255            elif nat == 'hydro':
    211256                md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
    … …  
    214259                md = checkfield(md, 'fieldname', 'materials.rho_freshwater', '>', 0)
    215260            else:
    216                 raise RuntimeError("materials checkconsistency error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
     261                raise RuntimeError('materials checkconsistency error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')')
    217262
    218263        return md
    … …  
    222267        #1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum
    223268        WriteData(fid, prefix, 'name', 'md.materials.nature', 'data', naturetointeger(self.nature), 'format', 'IntMat', 'mattype', 3)
    224         WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER, this can evolve if you have classes
    225 
    226         for i in range(len(self.nature)):
    227             nat = self.nature[i]
    228             if nat == 'ice':
    229                 WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 3, 'format', 'Integer')
     269        WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER: this can evolve if you have classes
     270        for i in range(len(self.nature)):
     271            nat = self.nature[i]
     272            if nat == 'ice':
    230273                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double')
    231274                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_water', 'format', 'Double')
    … …  
    236279                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'thermalconductivity', 'format', 'Double')
    237280                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'temperateiceconductivity', 'format', 'Double')
     281                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'effectiveconductivity_averaging', 'format', 'Integer')
    238282                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'meltingpoint', 'format', 'Double')
    239283                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'beta', 'format', 'Double')
    … …  
    251295                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'density', 'format', 'DoubleMat', 'mattype', 3)
    252296                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'viscosity', 'format', 'DoubleMat', 'mattype', 3)
    253                 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'isburgers', 'format', 'DoubleMat', 'mattype', 3)
     297                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheologymodel', 'format', 'DoubleMat', 'mattype', 3)
    254298                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_viscosity', 'format', 'DoubleMat', 'mattype', 3)
    255299                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_mu', 'format', 'DoubleMat', 'mattype', 3)
     300                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'ebm_alpha', 'format', 'DoubleMat', 'mattype', 3)
     301                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'ebm_delta', 'format', 'DoubleMat', 'mattype', 3)
     302                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'ebm_taul', 'format', 'DoubleMat', 'mattype', 3)
     303                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'ebm_tauh', 'format', 'DoubleMat', 'mattype', 3)
     304                # Compute earth density compatible with our layer density distribution
     305                earth_density = 0
     306                for i in range(self.numlayers):
     307                    earth_density = earth_density + (pow(self.radius[i + 1], 3) - pow(self.radius[i], 3)) * self.density[i]
     308                earth_density = earth_density / pow(self.radius[self.numlayers], 3)
     309                self.earth_density = earth_density
    256310            elif nat == 'hydro':
    257311                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double')
    258312                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_water', 'format', 'Double')
    259                 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
    260313                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_freshwater', 'format', 'Double')
    261314            else:
    262                 raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
     315                raise RuntimeError('materials constructor error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')')
     316        WriteData(fid, prefix, 'data', self.earth_density, 'name', 'md.materials.earth_density', 'format', 'Double')
    263317    #}}}
    264318
    … …  
    277331
    278332    for i in range(len(intnat)):
    279         if strnat[i] == 'damageice':
     333        str_nat = strnat[i]
     334        if str_nat == 'damageice':
    280335            intnat[i] = 1
    281         elif strnat[i] == 'estar':
     336        elif str_nat == 'estar':
    282337            intnat[i] = 2
    283         elif strnat[i] == 'ice':
     338        elif str_nat == 'ice':
    284339            intnat[i] = 3
    285         elif strnat[i] == 'enhancedice':
     340        elif str_nat == 'enhancedice':
    286341            intnat[i] = 4
    287         # elif strnat[i] == 'materials':
     342        # elif str_nat == 'materials':
    288343        #     intnat[i] = 5 #this case will never happen, kept to ensure equivalent of codes between IoCodeToMeterialsEnum and IoCodeToNatureEnum
    289         elif strnat[i] == 'litho':
     344        elif str_nat == 'litho':
    290345            intnat[i] = 6
    291         elif strnat[i] == 'hydro':
     346        elif str_nat == 'hydro':
    292347            intnat[i] = 7
    293348        else:
    294             raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
     349            raise RuntimeError('materials constructor error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')')
    295350
    296351    return intnat
Note: See TracChangeset for help on using the changeset viewer.