Changeset 25169


Ignore:
Timestamp:
06/27/20 01:49:04 (5 years ago)
Author:
jdquinn
Message:

CHG: More MATLAB -> Python translations; cleanup

Location:
issm/trunk-jpl
Files:
2 added
17 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/classes/dsl.py

    r24621 r25169  
    11import numpy as np
     2
     3from checkfield import checkfield
    24from fielddisplay import fielddisplay
    3 from checkfield import checkfield
     5from project3d import project3d
    46from WriteData import WriteData
    5 from project3d import project3d
    67
    78
    89class dsl(object):
    9     """
    10     dsl Class definition
     10    '''
     11    DSL - slass definition
    1112
    12        Usage:
    13           dsl = dsl()
    14     """
     13        Usage:
     14            dsl = dsl()
     15    '''
    1516
    16     def __init__(self):  # {{{
     17    def __init__(self): #{{{
    1718        self.global_average_thermosteric_sea_level_change = 0 #corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)
    1819        self.sea_surface_height_change_above_geoid = float('NaN') #corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)
    … …  
    2122    #}}}
    2223
    23     def __repr__(self):  # {{{
    24         string = "   dsl parameters:"
    25         string = "%s\n%s" % (string, fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)'))
    26         string = "%s\n%s" % (string, fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)'))
    27         string = "%s\n%s" % (string, fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)'))
    28         string = "%s\n%s" % (string, fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid'))
    29         return string
     24    def __repr__(self): #{{{
     25        s = "   dsl parameters:"
     26        s = "%s\n%s" % (s, fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)'))
     27        s = "%s\n%s" % (s, fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)'))
     28        s = "%s\n%s" % (s, fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)'))
     29        s = "%s\n%s" % (s, fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid'))
     30
     31        return s
    3032    #}}}
    3133
    32     def extrude(self, md):  # {{{
     34    def extrude(self, md): #{{{
    3335        self.sea_surface_height_change_above_geoid = project3d(md, 'vector', self.sea_surface_height_change_above_geoid, 'type', 'node')
    3436        self.sea_water_pressure_change_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_change_at_sea_floor, 'type', 'node')
    … …  
    3638    #}}}
    3739
    38     def defaultoutputs(self, md):  # {{{
     40    def defaultoutputs(self, md): #{{{
    3941        return []
    4042    #}}}
    4143
    42     def checkconsistency(self, md, solution, analyses):  # {{{
     44    def checkconsistency(self, md, solution, analyses): #{{{
    4345        # Early retun
    44         if not 'SealevelriseAnalysis' in analyses:
     46        if 'SealevelriseAnalysis' not in analyses:
    4547            return md
    4648
    … …  
    6062    # }}}
    6163
    62     def marshall(self, prefix, md, fid):    # {{{
     64    def marshall(self, prefix, md, fid):   #{{{
    6365        yts = md.constants.yts
    6466        WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer')

  • issm/trunk-jpl/src/m/classes/dslmme.m

    r24505 r25169  
    1515        end
    1616        methods
    17                 function self = extrude(self,md) % {{{
    18                         for i=1:length(self.global_average_thermosteric_sea_level_change),
    19                                 self.sea_surface_height_change_above_geoid{i}=project3d(md,'vector',self.sea_surface_height_change_above_geoid{i},'type','node','layer',1);
    20                                 self.sea_water_pressure_change_at_sea_floor{i}=project3d(md,'vector',self.sea_water_pressure_change_at_sea_floor{i},'type','node','layer',1);
    21                         end
    22                 end % }}}
    2317                function self = dslmme(varargin) % {{{
    2418                        switch nargin
    … …  
    8074
    8175                end % }}}
     76                function self = extrude(self,md) % {{{
     77                        for i=1:length(self.global_average_thermosteric_sea_level_change),
     78                                self.sea_surface_height_change_above_geoid{i}=project3d(md,'vector',self.sea_surface_height_change_above_geoid{i},'type','node','layer',1);
     79                                self.sea_water_pressure_change_at_sea_floor{i}=project3d(md,'vector',self.sea_water_pressure_change_at_sea_floor{i},'type','node','layer',1);
     80                        end
     81                end % }}}
    8282                function savemodeljs(self,fid,modelname) % {{{
    8383                       

  • issm/trunk-jpl/src/m/classes/geometry.m

    r24861 r25169  
    6161                                return;
    6262                        else
    63                                 md = checkfield(md,'fieldname','geometry.surface'  ,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
     63                                md = checkfield(md,'fieldname','geometry.surface' ,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    6464                                md = checkfield(md,'fieldname','geometry.base'      ,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    6565                                md = checkfield(md,'fieldname','geometry.thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);

  • issm/trunk-jpl/src/m/classes/geometry.py

    r24861 r25169  
    11import numpy as np
     2
     3from checkfield import checkfield
     4from fielddisplay import fielddisplay
    25from project3d import project3d
    3 from fielddisplay import fielddisplay
    4 from checkfield import checkfield
    56from WriteData import WriteData
    67
    78
    89class geometry(object):
    9     """
     10    '''
    1011    GEOMETRY class definition
    1112
    12        Usage:
    13           geometry = geometry()
    14     """
     13        Usage:
     14            geometry = geometry()
     15    '''
    1516
    16     def __init__(self):  # {{{
    17         self.surface = float('NaN')
    18         self.thickness = float('NaN')
    19         self.base = float('NaN')
    20         self.bed = float('NaN')
    21         self.hydrostatic_ratio = float('NaN')
     17    def __init__(self): #{{{
     18        self.surface = np.nan
     19        self.thickness = np.nan
     20        self.base = np.nan
     21        self.bed = np.nan
     22        self.hydrostatic_ratio = np.nan
    2223
    23     #set defaults
     24        #set defaults
    2425        self.setdefaultparameters()
    25 
    2626    #}}}
    2727
    28     def __repr__(self):  # {{{
    29 
    30         string = "   geometry parameters:"
    31         string = "%s\n%s" % (string, fielddisplay(self, 'surface', 'ice upper surface elevation [m]'))
    32         string = "%s\n%s" % (string, fielddisplay(self, 'thickness', 'ice thickness [m]'))
    33         string = "%s\n%s" % (string, fielddisplay(self, 'base', 'ice base elevation [m]'))
    34         string = "%s\n%s" % (string, fielddisplay(self, 'bed', 'bed elevation [m]'))
    35         return string
     28    def __repr__(self): #{{{
     29        s = "   geometry parameters:"
     30        s = "%s\n%s" % (s, fielddisplay(self, 'surface', 'ice upper surface elevation [m]'))
     31        s = "%s\n%s" % (s, fielddisplay(self, 'thickness', 'ice thickness [m]'))
     32        s = "%s\n%s" % (s, fielddisplay(self, 'base', 'ice base elevation [m]'))
     33        s = "%s\n%s" % (s, fielddisplay(self, 'bed', 'bed elevation [m]'))
     34       
     35        return s
    3636    #}}}
    3737
    38     def extrude(self, md):  # {{{
     38    def extrude(self, md): #{{{
    3939        self.surface = project3d(md, 'vector', self.surface, 'type', 'node')
    4040        self.thickness = project3d(md, 'vector', self.thickness, 'type', 'node')
    … …  
    4545    #}}}
    4646
    47     def setdefaultparameters(self):  # {{{
     47    def setdefaultparameters(self): #{{{
    4848        return self
    4949    #}}}
    5050
    51     def checkconsistency(self, md, solution, analyses):  # {{{
     51    def checkconsistency(self, md, solution, analyses): #{{{
    5252        if (solution == 'TransientSolution' and md.transient.isgia) or (solution == 'GiaSolution'):
    53             md = checkfield(md, 'fieldname', 'geometry.thickness', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
     53            md = checkfield(md, 'fieldname', 'geometry.thickness', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
     54        elif solution == 'SealevelriseSolution':
     55            md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    5456        elif solution == 'LoveSolution':
    5557            return
    … …  
    7274    # }}}
    7375
    74     def marshall(self, prefix, md, fid):  # {{{
     76    def marshall(self, prefix, md, fid): #{{{
    7577        WriteData(fid, prefix, 'object', self, 'fieldname', 'surface', 'format', 'DoubleMat', 'mattype', 1)
    7678        WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)

  • issm/trunk-jpl/src/m/classes/giamme.py

    r25165 r25169  
    5656
    5757        if isinstance(self.Ngia, list) or self.Ngia.ndim == 1: #list or 1D numpy.ndarray
    58             WriteData(fid, prefix, 'data', np.zeros(md.mesh.numberofvertices), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ngia')
    59             WriteData(fid, prefix, 'data', np.zeros(md.mesh.numberofvertices), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ugia')
     58            WriteData(fid, prefix, 'data', np.zeros((md.mesh.numberofvertices, 1)), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ngia')
     59            WriteData(fid, prefix, 'data', np.zeros((md.mesh.numberofvertices, 1)), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ugia')
    6060            WriteData(fid, prefix, 'data', 1, 'name', 'md.gia.modelid', 'format', 'Double')
    6161            WriteData(fid, prefix, 'name', 'md.gia.nummodels', 'data', 1, 'format', 'Integer')

  • issm/trunk-jpl/src/m/classes/matdamageice.m

    r23702 r25169  
    120120                        md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval'});
    121121                        md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]);
    122            
     122                       
    123123                        if ismember('GiaAnalysis',analyses),
    124124                                md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1);

  • issm/trunk-jpl/src/m/classes/matdamageice.py

    r25125 r25169  
    130130        md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
    131131        md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
     132
     133        if 'GiaAnalysis' in analyses:
     134            md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1)
     135            md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1)
     136            md = checkfield(md,'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1)
     137            md = checkfield(md,'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1)
     138
     139        if 'SealevelriseAnalysis' in analyses:
     140                md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', 1)
     141
    132142        return md
    133143    # }}}

  • issm/trunk-jpl/src/m/classes/matice.m

    r24741 r25169  
    6868
    6969                        %water viscosity (N.s/m^2)
    70                         self.mu_water=0.001787; 
     70                        self.mu_water=0.001787;
    7171
    7272                        %ice heat capacity cp (J/kg/K)
    … …  
    178178                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2);
    179179                        WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String');
    180 
    181180                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_shear_modulus','format','Double');
    182181                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_density','format','Double','scale',10^3);

  • issm/trunk-jpl/src/m/classes/matice.py

    r25125 r25169  
     1import numpy as np
     2
    13from fielddisplay import fielddisplay
    24from project3d import project3d
    … …  
    68
    79class matice(object):
    8     """
     10    '''
    911    MATICE class definition
    1012
    11        Usage:
    12           matice = matice()
    13     """
     13    Usage:
     14            matice = matice()
     15    '''
    1416
    15     def __init__(self):  # {{{
     17    def __init__(self): #{{{
    1618        self.rho_ice = 0.
    1719        self.rho_water = 0.
    … …  
    2729        self.mixed_layer_capacity = 0.
    2830        self.thermal_exchange_velocity = 0.
    29         self.rheology_B = float('NaN')
    30         self.rheology_n = float('NaN')
     31        self.rheology_B = np.nan
     32        self.rheology_n = np.nan
    3133        self.rheology_law = ''
    3234
    33         #giaivins:
     35        #giaivins
    3436        self.lithosphere_shear_modulus = 0.
    3537        self.lithosphere_density = 0.
    … …  
    3739        self.mantle_density = 0.
    3840
    39         #SLR
    40         self.earth_density = 5512
     41        #slr
     42        self.earth_density = 0
    4143        self.setdefaultparameters()
    4244    #}}}
    4345
    44     def __repr__(self):  # {{{
    45         string = "   Materials:"
     46    def __repr__(self): #{{{
     47        s = "   Materials:"
     48        s = "%s\n%s" % (s, fielddisplay(self, "rho_ice", "ice density [kg/m^3]"))
     49        s = "%s\n%s" % (s, fielddisplay(self, "rho_water", "water density [kg/m^3]"))
     50        s = "%s\n%s" % (s, fielddisplay(self, "rho_freshwater", "fresh water density [kg/m^3]"))
     51        s = "%s\n%s" % (s, fielddisplay(self, "mu_water", "water viscosity [Ns/m^2]"))
     52        s = "%s\n%s" % (s, fielddisplay(self, "heatcapacity", "heat capacity [J/kg/K]"))
     53        s = "%s\n%s" % (s, fielddisplay(self, "thermalconductivity", "ice thermal conductivity [W/m/K]"))
     54        s = "%s\n%s" % (s, fielddisplay(self, "temperateiceconductivity", "temperate ice thermal conductivity [W/m/K]"))
     55        s = "%s\n%s" % (s, fielddisplay(self, "effectiveconductivity_averaging", "computation of effectiveconductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)"))
     56        s = "%s\n%s" % (s, fielddisplay(self, "meltingpoint", "melting point of ice at 1atm in K"))
     57        s = "%s\n%s" % (s, fielddisplay(self, "latentheat", "latent heat of fusion [J/m^3]"))
     58        s = "%s\n%s" % (s, fielddisplay(self, "beta", "rate of change of melting point with pressure [K/Pa]"))
     59        s = "%s\n%s" % (s, fielddisplay(self, "mixed_layer_capacity", "mixed layer capacity [W/kg/K]"))
     60        s = "%s\n%s" % (s, fielddisplay(self, "thermal_exchange_velocity", "thermal exchange velocity [m/s]"))
     61        s = "%s\n%s" % (s, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1/n)]"))
     62        s = "%s\n%s" % (s, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
     63        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'"))
     64        s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
     65        s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_density", "Lithosphere density [g/cm^-3]"))
     66        s = "%s\n%s" % (s, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
     67        s = "%s\n%s" % (s, fielddisplay(self, "mantle_density", "Mantle density [g/cm^-3]"))
     68        s = "%s\n%s" % (s, fielddisplay(self, "earth_density", "Mantle density [kg/m^-3]"))
    4669
    47         string = "%s\n%s" % (string, fielddisplay(self, "rho_ice", "ice density [kg / m^3]"))
    48         string = "%s\n%s" % (string, fielddisplay(self, "rho_water", "water density [kg / m^3]"))
    49         string = "%s\n%s" % (string, fielddisplay(self, "rho_freshwater", "fresh water density [kg / m^3]"))
    50         string = "%s\n%s" % (string, fielddisplay(self, "mu_water", "water viscosity [N s / m^2]"))
    51         string = "%s\n%s" % (string, fielddisplay(self, "heatcapacity", "heat capacity [J / kg / K]"))
    52         string = "%s\n%s" % (string, fielddisplay(self, "thermalconductivity", "ice thermal conductivity [W / m / K]"))
    53         string = "%s\n%s" % (string, fielddisplay(self, "temperateiceconductivity", "temperate ice thermal conductivity [W / m / K]"))
    54         string = "%s\n%s" % (string, fielddisplay(self, "effectiveconductivity_averaging", "computation of effectiveconductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)"))
    55         string = "%s\n%s" % (string, fielddisplay(self, "meltingpoint", "melting point of ice at 1atm in K"))
    56         string = "%s\n%s" % (string, fielddisplay(self, "latentheat", "latent heat of fusion [J / m^3]"))
    57         string = "%s\n%s" % (string, fielddisplay(self, "beta", "rate of change of melting point with pressure [K / Pa]"))
    58         string = "%s\n%s" % (string, fielddisplay(self, "mixed_layer_capacity", "mixed layer capacity [W / kg / K]"))
    59         string = "%s\n%s" % (string, fielddisplay(self, "thermal_exchange_velocity", "thermal exchange velocity [m / s]"))
    60         string = "%s\n%s" % (string, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1 / n)]"))
    61         string = "%s\n%s" % (string, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
    62         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'"))
    63         string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
    64         string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_density", "Lithosphere density [g / cm^ - 3]"))
    65         string = "%s\n%s" % (string, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
    66         string = "%s\n%s" % (string, fielddisplay(self, "mantle_density", "Mantle density [g / cm^ - 3]"))
    67         string = "%s\n%s" % (string, fielddisplay(self, "earth_density", "Mantle density [kg / m^ - 3]"))
    68         return string
     70        return s
    6971    #}}}
    7072
    71     def extrude(self, md):  # {{{
     73    def extrude(self, md): #{{{
    7274        self.rheology_B = project3d(md, 'vector', self.rheology_B, 'type', 'node')
    7375        self.rheology_n = project3d(md, 'vector', self.rheology_n, 'type', 'element')
    … …  
    7577    #}}}
    7678
    77     def setdefaultparameters(self):  # {{{
    78         #ice density (kg / m^3)
     79    def setdefaultparameters(self): #{{{
     80        #ice density (kg/m^3)
    7981        self.rho_ice = 917.
    80         #ocean water density (kg / m^3)
     82        #ocean water density (kg/m^3)
    8183        self.rho_water = 1023.
    82         #fresh water density (kg / m^3)
     84        #fresh water density (kg/m^3)
    8385        self.rho_freshwater = 1000.
    84         #water viscosity (N.s / m^2)
     86        #water viscosity (N.s/m^2)
    8587        self.mu_water = 0.001787
    86         #ice heat capacity cp (J / kg / K)
     88        #ice heat capacity cp (J/kg/K)
    8789        self.heatcapacity = 2093.
    88         #ice latent heat of fusion L (J / kg)
     90        #ice latent heat of fusion L (J/kg)
    8991        self.latentheat = 3.34e5
    90         #ice thermal conductivity (W / m / K)
     92        #ice thermal conductivity (W/m/K)
    9193        self.thermalconductivity = 2.4
     94        #temperate ice thermal conductivity (W/m/K)
     95        self.temperateiceconductivity = 0.24
    9296        #computation of effective conductivity
    9397        self.effectiveconductivity_averaging = 1
    94         #temperate ice thermal conductivity (W / m / K)
    95         self.temperateiceconductivity = 0.24
    9698        #the melting point of ice at 1 atmosphere of pressure in K
    9799        self.meltingpoint = 273.15
    98         #rate of change of melting point with pressure (K / Pa)
     100        #rate of change of melting point with pressure (K/Pa)
    99101        self.beta = 9.8e-8
    100         #mixed layer (ice-water interface) heat capacity (J / kg / K)
     102        #mixed layer (ice-water interface) heat capacity (J/kg/K)
    101103        self.mixed_layer_capacity = 3974.
    102         #thermal exchange velocity (ice-water interface) (m / s)
     104        #thermal exchange velocity (ice-water interface) (m/s)
    103105        self.thermal_exchange_velocity = 1.00e-4
    104106        #Rheology law: what is the temperature dependence of B with T
    … …  
    107109
    108110        # GIA:
    109         self.lithosphere_shear_modulus = 6.7e10  # (Pa)
    110         self.lithosphere_density = 3.32  # (g / cm^ - 3)
     111        self.lithosphere_shear_modulus = 6.7e10 # (Pa)
     112        self.lithosphere_density = 3.32  # (g/cm^-3)
    111113        self.mantle_shear_modulus = 1.45e11  # (Pa)
    112         self.mantle_density = 3.34  # (g / cm^ - 3)
     114        self.mantle_density = 3.34  # (g/cm^-3)
    113115
    114         #SLR
    115         self.earth_density = 5512  # average density of the Earth, (kg / m^3)
    116         return self
     116        # SLR
     117        self.earth_density = 5512  # average density of the Earth, (kg/m^3)
    117118    #}}}
    118119
    119     def checkconsistency(self, md, solution, analyses):  # {{{
    120         md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
    121         md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0)
    122         md = checkfield(md, 'fieldname', 'materials.rho_freshwater', '>', 0)
    123         md = checkfield(md, 'fieldname', 'materials.mu_water', '>', 0)
    124         md = checkfield(md, 'fieldname', 'materials.rheology_B', '>', 0, 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    125         md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements])
    126         md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O'])
    127         md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
    128         md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', [1])
    129         md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', [1])
    130         md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', [1])
    131         md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
    132         md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
     120    def checkconsistency(self, md, solution, analyses): #{{{
     121        if solution != 'SealevelriseSolution':
     122            md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
     123            md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0)
     124            md = checkfield(md, 'fieldname', 'materials.rho_freshwater', '>', 0)
     125            md = checkfield(md, 'fieldname', 'materials.mu_water', '>', 0)
     126            md = checkfield(md, 'fieldname', 'materials.rheology_B', '>', 0, 'timeseries', 1, 'NaN', 1, 'Inf', 1)
     127            md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements])
     128            md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O'])
     129            md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
     130
     131        if 'GiaAnalysis' in analyses:
     132            md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', [1])
     133            md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', [1])
     134            md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', [1])
     135            md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
     136
     137        if 'SealevelriseAnalysis' in analyses:
     138            md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
     139
    133140        return md
    134     # }}}
     141    #}}}
    135142
    136     def marshall(self, prefix, md, fid):  # {{{
     143    def marshall(self, prefix, md, fid): #{{{
    137144        WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 3, 'format', 'Integer')
    138145        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double')
    … …  
    151158        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_B', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    152159        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2)
    153         WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
    154 
     160        WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 's')
    155161        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
    156162        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10.**3.)
    … …  
    158164        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10.**3.)
    159165        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
    160 
    161     # }}}
     166    #}}}

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

    r25168 r25169  
    5050
    5151    def __repr__(self):  # {{{
    52         string = '   slr parameters:'
    53         string = "%s\n%s" % (string, fielddisplay(self, 'deltathickness', 'thickness change: ice height equivalent [m]'))
    54         string = "%s\n%s" % (string, fielddisplay(self, 'sealevel', 'current sea level (prior to computation) [m]'))
    55         string = "%s\n%s" % (string, fielddisplay(self, 'spcthickness', 'thickness constraints (NaN means no constraint) [m]'))
    56         string = "%s\n%s" % (string, fielddisplay(self, 'reltol', 'sea level rise relative convergence criterion, (NaN: not applied)'))
    57         string = "%s\n%s" % (string, fielddisplay(self, 'abstol', 'sea level rise absolute convergence criterion, (default, NaN: not applied)'))
    58         string = "%s\n%s" % (string, fielddisplay(self, 'maxiter', 'maximum number of nonlinear iterations'))
    59         string = "%s\n%s" % (string, fielddisplay(self, 'love_h', 'load Love number for radial displacement'))
    60         string = "%s\n%s" % (string, fielddisplay(self, 'love_k', 'load Love number for gravitational potential perturbation'))
    61         string = "%s\n%s" % (string, fielddisplay(self, 'love_l', 'load Love number for horizontal displaements'))
    62         string = "%s\n%s" % (string, fielddisplay(self, 'tide_love_k', 'tidal load Love number (degree 2)'))
    63         string = "%s\n%s" % (string, fielddisplay(self, 'tide_love_h', 'tidal load Love number (degree 2)'))
    64         string = "%s\n%s" % (string, fielddisplay(self, 'fluid_love', 'secular fluid Love number'))
    65         string = "%s\n%s" % (string, fielddisplay(self, 'equatorial_moi', 'mean equatorial moment of inertia [kg m^2]'))
    66         string = "%s\n%s" % (string, fielddisplay(self, 'polar_moi', 'polar moment of inertia [kg m^2]'))
    67         string = "%s\n%s" % (string, fielddisplay(self, 'angular_velocity', 'mean rotational velocity of earth [per second]'))
    68         string = "%s\n%s" % (string, fielddisplay(self, 'ocean_area_scaling', 'correction for model representation of ocean area [default: No correction]'))
    69         string = "%s\n%s" % (string, fielddisplay(self, 'hydro_rate', 'rate of hydrological expansion [mm / yr]'))
    70         string = "%s\n%s" % (string, fielddisplay(self, 'geodetic', 'compute geodetic SLR? (in addition to steric?) default 0'))
    71         string = "%s\n%s" % (string, fielddisplay(self, 'geodetic_run_frequency', 'how many time steps we skip before we run SLR solver during transient (default: 1)'))
    72         string = "%s\n%s" % (string, fielddisplay(self, 'rigid', 'rigid earth graviational potential perturbation'))
    73         string = "%s\n%s" % (string, fielddisplay(self, 'elastic', 'elastic earth graviational potential perturbation'))
    74         string = "%s\n%s" % (string, fielddisplay(self, 'rotation', 'earth rotational potential perturbation'))
    75         string = "%s\n%s" % (string, fielddisplay(self, 'degacc', 'accuracy (default .01 deg) for numerical discretization of the Green''s functions'))
    76         string = "%s\n%s" % (string, fielddisplay(self, 'transitions', 'indices into parts of the mesh that will be icecaps'))
    77         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    78 
    79         return string
     52        s = '   slr parameters:'
     53        s = "%s\n%s" % (s, fielddisplay(self, 'deltathickness', 'thickness change: ice height equivalent [m]'))
     54        s = "%s\n%s" % (s, fielddisplay(self, 'sealevel', 'current sea level (prior to computation) [m]'))
     55        s = "%s\n%s" % (s, fielddisplay(self, 'spcthickness', 'thickness constraints (NaN means no constraint) [m]'))
     56        s = "%s\n%s" % (s, fielddisplay(self, 'reltol', 'sea level rise relative convergence criterion, (NaN: not applied)'))
     57        s = "%s\n%s" % (s, fielddisplay(self, 'abstol', 'sea level rise absolute convergence criterion, (default, NaN: not applied)'))
     58        s = "%s\n%s" % (s, fielddisplay(self, 'maxiter', 'maximum number of nonlinear iterations'))
     59        s = "%s\n%s" % (s, fielddisplay(self, 'love_h', 'load Love number for radial displacement'))
     60        s = "%s\n%s" % (s, fielddisplay(self, 'love_k', 'load Love number for gravitational potential perturbation'))
     61        s = "%s\n%s" % (s, fielddisplay(self, 'love_l', 'load Love number for horizontal displaements'))
     62        s = "%s\n%s" % (s, fielddisplay(self, 'tide_love_k', 'tidal load Love number (degree 2)'))
     63        s = "%s\n%s" % (s, fielddisplay(self, 'tide_love_h', 'tidal load Love number (degree 2)'))
     64        s = "%s\n%s" % (s, fielddisplay(self, 'fluid_love', 'secular fluid Love number'))
     65        s = "%s\n%s" % (s, fielddisplay(self, 'equatorial_moi', 'mean equatorial moment of inertia [kg m^2]'))
     66        s = "%s\n%s" % (s, fielddisplay(self, 'polar_moi', 'polar moment of inertia [kg m^2]'))
     67        s = "%s\n%s" % (s, fielddisplay(self, 'angular_velocity', 'mean rotational velocity of earth [per second]'))
     68        s = "%s\n%s" % (s, fielddisplay(self, 'ocean_area_scaling', 'correction for model representation of ocean area [default: No correction]'))
     69        s = "%s\n%s" % (s, fielddisplay(self, 'hydro_rate', 'rate of hydrological expansion [mm / yr]'))
     70        s = "%s\n%s" % (s, fielddisplay(self, 'geodetic', 'compute geodetic SLR? (in addition to steric?) default 0'))
     71        s = "%s\n%s" % (s, fielddisplay(self, 'geodetic_run_frequency', 'how many time steps we skip before we run SLR solver during transient (default: 1)'))
     72        s = "%s\n%s" % (s, fielddisplay(self, 'rigid', 'rigid earth graviational potential perturbation'))
     73        s = "%s\n%s" % (s, fielddisplay(self, 'elastic', 'elastic earth graviational potential perturbation'))
     74        s = "%s\n%s" % (s, fielddisplay(self, 'rotation', 'earth rotational potential perturbation'))
     75        s = "%s\n%s" % (s, fielddisplay(self, 'degacc', 'accuracy (default .01 deg) for numerical discretization of the Green''s functions'))
     76        s = "%s\n%s" % (s, fielddisplay(self, 'transitions', 'indices into parts of the mesh that will be icecaps'))
     77        s = "%s\n%s" % (s, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
     78
     79        return s
    8080    # }}}
    8181
    … …  
    143143        md = checkfield(md, 'fieldname', 'slr.hydro_rate', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    144144        md = checkfield(md, 'fieldname', 'slr.degacc', 'size', [1, 1], '>=', 1e-10)
    145         md = checkfield(md, 'fieldname', 'slr.requested_outputs', 'stringrow', 1)
     145        md = checkfield(md, 'fieldname', 'slr.requested_outputs', 'srow', 1)
    146146        md = checkfield(md, 'fieldname', 'slr.horiz', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
    147147
    … …  
    202202            outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
    203203            outputs = outputscopy
    204         WriteData(fid, prefix, 'data', outputs, 'name', 'md.slr.requested_outputs', 'format', 'StringArray')
    205     # }}}
     204        WriteData(fid, prefix, 'data', outputs, 'name', 'md.slr.requested_outputs', 'format', 'sArray')
     205    # }}}

  • issm/trunk-jpl/src/m/classes/solidearth.py

    r25161 r25169  
    2424        self.settings = solidearthsettings()
    2525        self.surfaceload = surfaceload()
    26         self.love = lovenumbers()
     26        self.lovenumbers = lovenumbers()
    2727        self.rotational = rotational()
    2828        self.planetradius = planetradius('earth')
    … …  
    6565        self.settings.checkconsistency(md, solution, analyses)
    6666        self.surfaceload.checkconsistency(md, solution, analyses)
    67         self.love.checkconsistency(md, solution, analyses)
     67        self.lovenumbers.checkconsistency(md, solution, analyses)
    6868        self.rotational.checkconsistency(md, solution, analyses)
    6969
    … …  
    8282        self.settings.marshall(prefix, md, fid)
    8383        self.surfaceload.marshall(prefix, md, fid)
    84         self.love.marshall(prefix, md, fid)
     84        self.lovenumbers.marshall(prefix, md, fid)
    8585        self.rotational.marshall(prefix, md, fid)
    8686

  • issm/trunk-jpl/src/m/classes/toolkits.m

    r24746 r25169  
    55
    66classdef toolkits < dynamicprops
    7     properties (SetAccess=public)
    8                  DefaultAnalysis  = struct();
    9                  RecoveryAnalysis = struct();
    10                  %The other properties are dynamic
    11          end
    12          methods (Static)
    13                  function self = loadobj(self) % {{{
    14                          % This function is directly called by matlab when a model object is
    15                          % loaded. Update old properties here
     7        properties (SetAccess=public)
     8                DefaultAnalysis  = struct();
     9                RecoveryAnalysis = struct();
     10                %The other properties are dynamic
     11        end
     12        methods (Static)
     13                function self = loadobj(self) % {{{
     14                        % This function is directly called by matlab when a model object is
     15                        % loaded. Update old properties here
    1616
    17                          if isempty(fieldnames(self.RecoveryAnalysis));
    18                                  disp('WARNING: updating toolkits (RecoveryAnalysis not set)');
    19                                  self.RecoveryAnalysis  = self.DefaultAnalysis;
    20                          end
    21                  end% }}}
    22          end
    23          methods
    24                  function self = toolkits(varargin) % {{{
    25                          switch nargin
    26                                  case 0
    27                                          self=setdefaultparameters(self);
    28                                  case 1
    29                                          self=structtoobj(self,varargin{1});
    30                                  otherwise
    31                                          error('constructor not supported');
    32                                  end
    33                          end % }}}
    34                  function self = addoptions(self,analysis,varargin) % {{{
    35                  % Usage example:
    36                  %    md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis',FSoptions());
    37                  %    md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis');
     17                        if isempty(fieldnames(self.RecoveryAnalysis));
     18                                disp('WARNING: updating toolkits (RecoveryAnalysis not set)');
     19                                self.RecoveryAnalysis  = self.DefaultAnalysis;
     20                        end
     21                end% }}}
     22        end
     23        methods
     24                function self = toolkits(varargin) % {{{
     25                        switch nargin
     26                                case 0
     27                                        self=setdefaultparameters(self);
     28                                case 1
     29                                        self=structtoobj(self,varargin{1});
     30                                otherwise
     31                                        error('constructor not supported');
     32                                end
     33                        end % }}}
     34                function self = addoptions(self,analysis,varargin) % {{{
     35                %ADDOPTIONS - add analysis to md.toolkits.analysis
     36                %
     37                %   Optional third parameter adds toolkits options to analysis.
     38                %
     39                %   Usage:
     40                %      md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis',FSoptions());
     41                %      md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis');
    3842
    39                          %Create dynamic property if property does not exist yet
    40                          if ~ismember(analysis,properties(self)),
    41                                  self.addprop(analysis);
    42                          end
     43                        %Create dynamic property if property does not exist yet
     44                        if ~ismember(analysis,properties(self)),
     45                                self.addprop(analysis);
     46                        end
    4347
    44                          %Add toolkits options to analysis
    45                          if nargin==3, self.(analysis) = varargin{1}; end
    46                  end
    47                  %}}}
    48                  function self = setdefaultparameters(self) % {{{
     48                        %Add toolkits options to analysis
     49                        if nargin==3,
     50                                self.(analysis) = varargin{1};
     51                        end
     52                end
     53                %}}}
     54                function self = setdefaultparameters(self) % {{{
    4955
    50                          %default toolkits:
    51                          if IssmConfig('_HAVE_PETSC_'),
    52                                  %MUMPS is the default toolkits
    53                                  if IssmConfig('_HAVE_MUMPS_'),
    54                                          self.DefaultAnalysis           = mumpsoptions();
    55                                  else
    56                                          self.DefaultAnalysis           = iluasmoptions();
    57                                  end
    58                          else
    59                                  if IssmConfig('_HAVE_MUMPS_'),
    60                                          self.DefaultAnalysis           = issmmumpssolver();
    61                                  elseif IssmConfig('_HAVE_GSL_'),
    62                                          self.DefaultAnalysis           = issmgslsolver();
    63                                  else
    64                                          disp('WARNING: Need at least Mumps or Gsl to define an issm solver type, no default solver assigned');
    65                                  end
    66                          end
     56                        %default toolkits:
     57                        if IssmConfig('_HAVE_PETSC_'),
     58                                %MUMPS is the default toolkits
     59                                if IssmConfig('_HAVE_MUMPS_'),
     60                                        self.DefaultAnalysis           = mumpsoptions();
     61                                else
     62                                        self.DefaultAnalysis           = iluasmoptions();
     63                                end
     64                        else
     65                                if IssmConfig('_HAVE_MUMPS_'),
     66                                        self.DefaultAnalysis           = issmmumpssolver();
     67                                elseif IssmConfig('_HAVE_GSL_'),
     68                                        self.DefaultAnalysis           = issmgslsolver();
     69                                else
     70                                        disp('WARNING: Need at least Mumps or Gsl to define an issm solver type, no default solver assigned');
     71                                end
     72                        end
    6773
    68                          %Use same solver for Recovery mode
    69                          self.RecoveryAnalysis = self.DefaultAnalysis;
     74                        %Use same solver for Recovery mode
     75                        self.RecoveryAnalysis = self.DefaultAnalysis;
    7076
    7177
    72                  end % }}}
    73                  function disp(self) % {{{
    74                          analyses=properties(self);
    75                          disp(sprintf('List of toolkits options per analysis:\n'));
    76                          for i=1:numel(analyses),
    77                                  analysis=analyses{i};
    78                                  disp([analysis ':']);
    79                                  disp(self.(analysis));
    80                          end
    81                  end % }}}
    82                  function md = checkconsistency(self,md,solution,analyses) % {{{
    83                          analyses=properties(self);
    84                          for i=1:numel(analyses),
    85                                  switch analyses{i}
    86                                          case 'DefaultAnalysis'
    87                                          case 'RecoveryAnalysis'
    88                                          case 'StressbalanceAnalysis'
    89                                          case 'GLheightadvectionAnalysis'
    90                                          case 'MasstransportAnalysis'
    91                                          case 'ThermalAnalysis'
    92                                          case 'EnthalpyAnalysis'
    93                                          case 'AdjointBalancethicknessAnalysis'
    94                                          case 'BalancethicknessAnalysis'
    95                                          case 'Balancethickness2Analysis'
    96                                          case 'BalancethicknessSoftAnalysis'
    97                                          case 'BalancevelocityAnalysis'
    98                                          case 'DamageEvolutionAnalysis'
    99                                          case 'LoveAnalysis'
    100                                          case 'EsaAnalysis'
    101                                          case 'SealevelriseAnalysis'
    102                                          otherwise
     78                end % }}}
     79                function disp(self) % {{{
     80                        analyses=properties(self);
     81                        disp(sprintf('List of toolkits options per analysis:\n'));
     82                        for i=1:numel(analyses),
     83                                analysis=analyses{i};
     84                                disp([analysis ':']);
     85                                disp(self.(analysis));
     86                        end
     87                end % }}}
     88                function md = checkconsistency(self,md,solution,analyses) % {{{
     89                        analyses=properties(self);
     90                        for i=1:numel(analyses),
     91                                switch analyses{i}
     92                                        case 'DefaultAnalysis'
     93                                        case 'RecoveryAnalysis'
     94                                        case 'StressbalanceAnalysis'
     95                                        case 'GLheightadvectionAnalysis'
     96                                        case 'MasstransportAnalysis'
     97                                        case 'ThermalAnalysis'
     98                                        case 'EnthalpyAnalysis'
     99                                        case 'AdjointBalancethicknessAnalysis'
     100                                        case 'BalancethicknessAnalysis'
     101                                        case 'Balancethickness2Analysis'
     102                                        case 'BalancethicknessSoftAnalysis'
     103                                        case 'BalancevelocityAnalysis'
     104                                        case 'DamageEvolutionAnalysis'
     105                                        case 'LoveAnalysis'
     106                                        case 'EsaAnalysis'
     107                                        case 'SealevelriseAnalysis'
     108                                        otherwise
    103109                                                md = checkmessage(md,['md.toolkits.' analyses{i} ' not supported yet']);
    104                                  end
    105                                  if isempty(fieldnames(self.(analyses{i})))
    106                                          md = checkmessage(md,['md.toolkits.' analyses{i} ' is empty']);
    107                                  end
    108                          end
    109                  end % }}}
    110                  function ToolkitsFile(toolkits,filename) % {{{
    111                          %TOOLKITSFILE - build toolkits file
    112                          %
    113                          %   Build a Petsc compatible options file, from the toolkits model field  + return options string.
    114                          %   This file will also be used when the toolkit used is 'issm' instead of 'petsc'
    115                          %
    116                          %   Usage:     ToolkitsFile(toolkits,filename);
     110                                end
     111                                if isempty(fieldnames(self.(analyses{i})))
     112                                        md = checkmessage(md,['md.toolkits.' analyses{i} ' is empty']);
     113                                end
     114                        end
     115                end % }}}
     116                function ToolkitsFile(toolkits,filename) % {{{
     117                        %TOOLKITSFILE - build toolkits file
     118                        %
     119                        %   Build a Petsc compatible options file, from the toolkits model field  + return options string.
     120                        %   This file will also be used when the toolkit used is 'issm' instead of 'petsc'
     121                        %
     122                        %   Usage:     ToolkitsFile(toolkits,filename);
    117123
    118                          %open file for writing
    119                          fid=fopen(filename,'w');
    120                          if fid==-1,
    121                                  error(['ToolkitsFile error: could not open ' filename ' for writing']);
    122                          end
     124                        %open file for writing
     125                        fid=fopen(filename,'w');
     126                        if fid==-1,
     127                                error(['ToolkitsFile error: could not open ' filename ' for writing']);
     128                        end
    123129
    124                          %write header
    125                          fprintf(fid,'%s%s%s\n','%Toolkits options file: ',filename,' written from Matlab toolkits array');
     130                        %write header
     131                        fprintf(fid,'%s%s%s\n','%Toolkits options file: ',filename,' written from Matlab toolkits array');
    126132
    127                          %start writing options
    128                          analyses=properties(toolkits);
    129                          for i=1:numel(analyses),
    130                                  analysis=analyses{i};
    131                                  options=toolkits.(analysis);
     133                        %start writing options
     134                        analyses=properties(toolkits);
     135                        for i=1:numel(analyses),
     136                                analysis=analyses{i};
     137                                options=toolkits.(analysis);
    132138
    133                                  %first write analysis:
    134                                  fprintf(fid,'\n+%s\n',analysis); %append a + to recognize it's an analysis enum
     139                                %first write analysis:
     140                                fprintf(fid,'\n+%s\n',analysis); %append a + to recognize it's an analysis enum
    135141
    136                                  %now, write options
    137                                  optionslist=fieldnames(options);
    138                                  for j=1:numel(optionslist),
    139                                          optionname=optionslist{j};
    140                                          optionvalue=options.(optionname);
     142                                %now, write options
     143                                optionslist=fieldnames(options);
     144                                for j=1:numel(optionslist),
     145                                        optionname=optionslist{j};
     146                                        optionvalue=options.(optionname);
    141147
    142                                          if isempty(optionvalue),
    143                                                  %this option has only one argument
    144                                                  fprintf(fid,'-%s\n',optionname);
    145                                          else
    146                                                  %option with value. value can be string or scalar
    147                                                  if isnumeric(optionvalue),
    148                                                          fprintf(fid,'-%s %g\n',optionname,optionvalue);
    149                                                  elseif ischar(optionvalue),
    150                                                          fprintf(fid,'-%s %s\n',optionname,optionvalue);
    151                                                  else
    152                                                          error(['ToolkitsFile error: option ' optionname ' is not well formatted']);
    153                                                  end
    154                                          end
    155                                  end
    156                          end
     148                                        if isempty(optionvalue),
     149                                                %this option has only one argument
     150                                                fprintf(fid,'-%s\n',optionname);
     151                                        else
     152                                                %option with value. value can be string or scalar
     153                                                if isnumeric(optionvalue),
     154                                                        fprintf(fid,'-%s %g\n',optionname,optionvalue);
     155                                                elseif ischar(optionvalue),
     156                                                        fprintf(fid,'-%s %s\n',optionname,optionvalue);
     157                                                else
     158                                                        error(['ToolkitsFile error: option ' optionname ' is not well formatted']);
     159                                                end
     160                                        end
     161                                end
     162                        end
    157163
    158                          fclose(fid);
    159                  end %}}}
    160                  function savemodeljs(self,fid,modelname) % {{{
     164                        fclose(fid);
     165                end %}}}
     166                function savemodeljs(self,fid,modelname) % {{{
    161167
    162                          analyses=properties(self);
    163                          for i=1:numel(analyses),
    164                                  if isempty(fieldnames(self.(analyses{i})))
    165                                          error(['md.toolkits.' analyses{i} ' is empty']);
    166                                  else
    167                                          writejsstruct(fid,[modelname '.toolkits.' analyses{i}],self.(analyses{i}));
    168                                  end
    169                          end
     168                        analyses=properties(self);
     169                        for i=1:numel(analyses),
     170                                if isempty(fieldnames(self.(analyses{i})))
     171                                        error(['md.toolkits.' analyses{i} ' is empty']);
     172                                else
     173                                        writejsstruct(fid,[modelname '.toolkits.' analyses{i}],self.(analyses{i}));
     174                                end
     175                        end
    170176
    171                  end % }}}
    172          end
     177                end % }}}
     178        end
    173179 end

  • issm/trunk-jpl/src/m/classes/toolkits.py

    r24213 r25169  
    88
    99class toolkits(object):
    10     """
     10    '''
    1111    TOOLKITS class definition
    1212
    13        Usage:
    14           self = toolkits()
    15     """
     13        Usage:
     14            self = toolkits()
     15    '''
    1616
    17     def __init__(self):  # {{{
     17    def __init__(self): #{{{
    1818        #default toolkits
    1919        if IssmConfig('_HAVE_PETSC_')[0]:
    … …  
    3434        self.RecoveryAnalysis = self.DefaultAnalysis
    3535
    36     #The other properties are dynamic
    37     # }}}
     36        #The other properties are dynamic
     37    #}}}
    3838
    39     def __repr__(self):  # {{{
     39    def __repr__(self): #{{{
    4040        s = "List of toolkits options per analysis:\n\n"
    4141        for analysis in list(vars(self).keys()):
    42             s += "%s\n" % fielddisplay(self, analysis, '')
     42            s += "{}\n".format(fielddisplay(self, analysis, ''))
    4343
    44             return s
    45     # }}}
     44        return s
     45    #}}}
    4646
    47     def addoptions(self, analysis, *args):  # {{{
     47    def addoptions(self, analysis, *args): #{{{
    4848        # Usage example:
    4949        #    md.toolkits = addoptions(md.toolkits, 'StressbalanceAnalysis', FSoptions())
    … …  
    5959
    6060        return self
    61     # }}}
     61    #}}}
    6262
    63     def checkconsistency(self, md, solution, analyses):  # {{{
     63    def checkconsistency(self, md, solution, analyses): #{{{
     64        # TODO
     65        # - Implement something closer to a switch as in
     66        # src/m/classes/toolkits.m?
     67        #
    6468        for analysis in list(vars(self).keys()):
    6569            if not getattr(self, analysis):
    66                 md.checkmessage("md.toolkits.%s is empty" % analysis)
     70                md.checkmessage("md.toolkits.{} is empty".format(analysis))
    6771
    6872        return md
    69     # }}}
     73    #}}}
    7074
    71     def ToolkitsFile(self, filename):  # {{{
    72         """
     75    def ToolkitsFile(self, filename): #{{{
     76        '''
    7377        TOOLKITSFILE - build toolkits file
    7478
    75            Build a Petsc compatible options file, from the toolkits model field + return options string
    76            This file will also be used when the toolkit used is 'issm' instead of 'petsc'
     79            Build a Petsc compatible options file, from the toolkits model
     80            field + return options string.
     81            This file will also be used when the toolkit used is 'issm' instead
     82            of 'petsc'.s
    7783
     84            Usage:
     85                ToolkitsFile(toolkits, filename)
     86        '''
    7887
    79            Usage:     ToolkitsFile(toolkits, filename)
    80         """
    81 
    82     #open file for writing
     88        #open file for writing
    8389        try:
    8490            fid = open(filename, 'w')
    … …  
    8692            raise IOError("ToolkitsFile error: could not open {}' for writing due to".format(filename), e)
    8793
    88     #write header
     94        #write header
    8995        fid.write("%s%s%s\n" % ('%Toolkits options file: ', filename, ' written from Python toolkits array'))
    9096
    91     #start writing options
     97        #start writing options
    9298        for analysis in list(vars(self).keys()):
    9399            options = getattr(self, analysis)
    94100
    95     #first write analysis:
     101            #first write analysis:
    96102            fid.write("\n+{}\n".format(analysis))  #append a + to recognize it's an analysis enum
    97     #now, write options
     103            #now, write options
    98104            for optionname, optionvalue in list(options.items()):
    99105
    … …  
    111117
    112118        fid.close()
    113     # }}}
     119    #}}}

  • issm/trunk-jpl/src/m/consistency/checkfield.py

    r24808 r25169  
    88
    99def checkfield(md, *args):
    10     """
     10    '''
    1111    CHECKFIELD - check field consistency
    1212
    13        Used to check model consistency.,
    14        Requires:
    15        'field' or 'fieldname' option. If 'fieldname' is provided, it will retrieve it from the model md. (md.(fieldname))
    16              If 'field' is provided, it will assume the argument following 'field' is a numeric array.
    17 
    18        Available options:
    19  - NaN: 1 if check that there is no NaN
    20  - size: [lines cols], NaN for non checked dimensions, or 'universal' for any input type (nodal, element, time series, etc)
    21  -> :  greater than provided value
    22  ->= : greater or equal to provided value
    23  - < :  smallerthan provided value
    24  - <=: smaller or equal to provided value
    25  - < vec:  smallerthan provided values on each vertex
    26  - timeseries: 1 if check time series consistency (size and time)
    27  - values: cell of strings or vector of acceptable values
    28  - numel: list of acceptable number of elements
    29  - cell: 1 if check that is cell
    30  - empty: 1 if check that non empty
    31  - message: overloaded error message
    32 
    33        Usage:
    34           md = checkfield(md, fieldname, options)
    35     """
     13        Used to check model consistency
     14       
     15        Requires:
     16            'field' or 'fieldname' option. If 'fieldname' is provided, it will retrieve it from the model md. (md.(fieldname))
     17            If 'field' is provided, it will assume the argument following
     18            'field' is a numeric array.
     19
     20        Available options:
     21        - NaN: 1 if check that there is no NaN
     22        - size: [lines cols], NaN for non checked dimensions, or 'universal'
     23        for any input type (nodal, element, time series, etc)
     24        - > :  greater than provided value
     25        - >= : greater or equal to provided value
     26        - < :  smallerthan provided value
     27        - <=: smaller or equal to provided value
     28        - < vec:  smallerthan provided values on each vertex
     29        - timeseries: 1 if check time series consistency (size and time)
     30        - values: list of acceptable values
     31        - numel: list of acceptable number of elements
     32        - cell: 1 if check that is cell
     33        - empty: 1 if check that non empty
     34        - message: overloaded error message
     35
     36        Usage:
     37            md = checkfield(md, fieldname, options)
     38    '''
    3639
    3740    #get options
    … …  
    143146    if options.getfieldvalue('cell', 0):
    144147        if not isinstance(field, (tuple, list, dict)):
    145             md = md.checkmessage(options.getfieldvalue('message', "field '{}' should be a cell".format(fieldname)))
     148            md = md.checkmessage(options.getfieldvalue('message', "field '{}' should be a tuple, list, or dict".format(fieldname)))
    146149
    147150    #check values

  • issm/trunk-jpl/src/m/solve/solve.py

    r25168 r25169  
    22import os
    33import shutil
     4
     5from ismodelselfconsistent import ismodelselfconsistent
     6from loadresultsfromcluster import loadresultsfromcluster
     7from marshall import marshall
    48from pairoptions import pairoptions
    5 from ismodelselfconsistent import ismodelselfconsistent
    6 from marshall import marshall
     9from preqmu import *
    710from waitonlock import waitonlock
    8 from loadresultsfromcluster import loadresultsfromcluster
    9 from preqmu import *
    10 #from MatlabFuncs import *
    1111
    1212
    … …  
    1515    SOLVE - apply solution sequence for this model
    1616
    17        Usage:
    18           md = solve(md, solutionstring, varargin)
    19           where varargin is a list of paired arguments of string OR enums
     17        Usage:
     18            md = solve(md, solutionstring, varargin)
     19            where varargin is a list of paired arguments of string OR enums
    2020
    2121        solution types available comprise:
    22          - 'Stressbalance'    or 'sb'
    23          - 'Masstransport'    or 'mt'
    24          - 'Thermal'          or 'th'
    25          - 'Steadystate'      or 'ss'
    26          - 'Transient'        or 'tr'
    27          - 'Balancethickness' or 'mc'
    28          - 'Balancevelocity'  or 'bv'
    29          - 'BedSlope'         or 'bsl'
    30          - 'SurfaceSlope'     or 'ssl'
    31          - 'Hydrology'        or 'hy'
    32          - 'DamageEvolution'  or 'da'
    33          - 'Gia'              or 'gia'
    34          - 'Esa'          or 'esa'
    35          - 'Sealevelrise'     or 'slr'
    36          - 'Love'             or 'lv'
     22            - 'Stressbalance'      or 'sb'
     23            - 'Masstransport'      or 'mt'
     24            - 'Thermal'            or 'th'
     25            - 'Steadystate'        or 'ss'
     26            - 'Transient'          or 'tr'
     27            - 'Balancethickness'  or 'mc'
     28            - 'Balancevelocity'    or 'bv'
     29            - 'BedSlope'           or 'bsl'
     30            - 'SurfaceSlope'       or 'ssl'
     31            - 'Hydrology'          or 'hy'
     32            - 'DamageEvolution'    or 'da'
     33            - 'Gia'                or 'gia'
     34            - 'Esa'                or 'esa'
     35            - 'Sealevelrise'       or 'slr'
     36            - 'Love'               or 'lv'
    3737
    38        extra options:
    39  - loadonly : does not solve. only load results
    40          - checkconsistency : 'yes' or 'no' (default is 'yes'), ensures checks on consistency of model
    41          - restart: 'directory name (relative to the execution directory) where the restart file is located.
     38        extra options:
     39            - loadonly          : does not solve. only load results
     40            - checkconsistency  : 'yes' or 'no' (default is 'yes'), ensures
     41                                  checks on consistency of model
     42            - restart           : 'directory name (relative to the execution
     43                                  directory) where the restart file is located
    4244
    43        Examples:
    44           md = solve(md, 'Stressbalance')
    45          md = solve(md, 'sb')
     45        Examples:
     46            md = solve(md, 'Stressbalance')
     47            md = solve(md, 'sb')
    4648    '''
    4749

  • issm/trunk-jpl/test/NightlyRun/test2002.m

    r25166 r25169  
    6565
    6666%eustatic run:
    67 md.solidearth.settings.rigid=0; md.solidearth.settings.elastic=0;md.solidearth.settings.rotation=0;
     67md.solidearth.settings.rigid=0;
     68md.solidearth.settings.elastic=0;
     69md.solidearth.settings.rotation=0;
    6870md=solve(md,'Sealevelrise');
    6971Seustatic=md.results.SealevelriseSolution.Sealevel;

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

    r25166 r25169  
    3939#mask:  {{{
    4040mask = gmtmask(md.mesh.lat, md.mesh.long)
    41 icemask = np.ones(md.mesh.numberofvertices)
     41icemask = np.ones((md.mesh.numberofvertices, 1))
    4242pos = np.where(mask == 0)[0]
    4343icemask[pos] = -1
    … …  
    4848
    4949#make sure that the ice level set is all inclusive:
    50 md.mask.land_levelset = np.zeros((md.mesh.numberofvertices))
    51 md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices))
     50md.mask.land_levelset = np.zeros((md.mesh.numberofvertices, 1))
     51md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, 1))
    5252
    5353#make sure that the elements that have loads are fully grounded
Note: See TracChangeset for help on using the changeset viewer.