Changeset 26301


Ignore:
Timestamp:
06/07/21 20:07:14 (4 years ago)
Author:
jdquinn
Message:

CHG: test2002 MATLAB -> Python; clean up

Location:
issm/trunk-jpl
Files:
1 added
29 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py

    r25253 r26301  
    55
    66def getlovenumbers(*args): #{{{
    7     '''
    8     GETLOVENUMBERS - provide love numbers retrieved from:
     7    """GETLOVENUMBERS - provide love numbers retrieved from:
    98    http://www.srosat.com/iag-jsg/loveNb.php in a chosen reference frame
    109
    11         Usage:
    12             series = love_numbers('type', 'loadingverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', 1001)
     10    Usage:
     11        series = love_numbers('type', 'loadingverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', 1001)
    1312
    14         - type = one of 'loadingverticaldisplacement',
    15         'loadinggravitationalpotential', 'loadinghorizontaldisplacement',
    16         'tidalverticaldisplacement', 'tidalgravitationalpotential',
    17         'tidalhorizontaldisplacement'
    18         - reference_frame = one of 'CM' (default) and 'CF'
    19         - maxdeg = default 1000
     13    - type = one of 'loadingverticaldisplacement',
     14    'loadinggravitationalpotential', 'loadinghorizontaldisplacement',
     15    'tidalverticaldisplacement', 'tidalgravitationalpotential',
     16    'tidalhorizontaldisplacement'
     17    - reference_frame = one of 'CM' (default) and 'CF'
     18    - maxdeg = default 1000
    2019
    21         Example:
    22             h = love_number
    23     '''
     20    Example:
     21        h = getlovenumbers('type', 'loadingverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg)
     22        k = getlovenumbers('type', 'loadinggravitationalpotential', 'referenceframe', 'CM', 'maxdeg', maxdeg)
     23        l = getlovenumbers('type', 'loadinghorizontaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg)
     24        th = getlovenumbers('type', 'tidalverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg)
     25        tk = getlovenumbers('type', 'tidalgravitationalpotential', 'referenceframe', 'CM', 'maxdeg', maxdeg)
     26        tl = getlovenumbers('type', 'tidalhorizontaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg)
     27    """
    2428
    2529    TYPES = [
    … …  
    3236    ]
    3337
    34     #recover options
     38    # Recover options
    3539    options = pairoptions(*args)
    3640    type = options.getfieldvalue('type')
    3741    frame = options.getfieldvalue('referenceframe', 'CM')
    3842    maxdeg = options.getfieldvalue('maxdeg', 1000)
     43
     44    if (maxdeg > 10000):
     45        raise Exception('PREM love numbers computed only for deg < 10,000. Request lower maxdeg')
    3946
    4047    if type not in TYPES:
    … …  
    1004510052    ])
    1004610053
    10047     #cut off series at maxdeg
     10054    # Cut off series at maxdeg
    1004810055    love_numbers = np.delete(love_numbers, range(maxdeg + 1, len(love_numbers)), axis=0)
    1004910056
    10050     #retrieve right type
     10057    # Retrieve right type
    1005110058    if type == 'loadingverticaldisplacement':
    1005210059        series = love_numbers[:, 0]
    … …  
    1006210069        series = love_numbers[:, 5]
    1006310070    else:
    10064         raise Exception("love_numbers error message: unknown type: {}".format(type))
     10071        raise Exception('love_numbers error message: unknown type: {}'.format(type))
    1006510072
    1006610073    #choose degree 1 term for CF reference system
    … …  
    1007510082            series[1] = 0.134
    1007610083    else:
    10077         raise Exception("love_numbers error message: unknown reference frame: {}".format(frame))
     10084        raise Exception('love_numbers error message: unknown reference frame: {}'.format(frame))
    1007810085
    1007910086    return series

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

    r22885 r26301  
    8484                        self.gradation=1.5;
    8585
    86                         %other criterias
     86                        %other criteria
    8787                        self.groundingline_resolution=500.;
    8888                        self.groundingline_distance=0.;

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

    r24213 r26301  
    55
    66class amr(object):
    7     """
    8     AMR Class definition
     7    """AMR Class definition
    98
    109    Usage:
    … …  
    1211    """
    1312
    14     def __init__(self):  # {{{
    15         self.hmin = 0.
    16         self.hmax = 0.
     13    def __init__(self): #{{{
     14        self.hmin = 0
     15        self.hmax = 0
    1716        self.fieldname = ''
    18         self.err = 0.
    19         self.keepmetric = 0.
    20         self.gradation = 0.
    21         self.groundingline_resolution = 0.
    22         self.groundingline_distance = 0.
    23         self.icefront_resolution = 0.
    24         self.icefront_distance = 0.
    25         self.thicknesserror_resolution = 0.
    26         self.thicknesserror_threshold = 0.
    27         self.thicknesserror_groupthreshold = 0.
    28         self.thicknesserror_maximum = 0.
    29         self.deviatoricerror_resolution = 0.
    30         self.deviatoricerror_threshold = 0.
    31         self.deviatoricerror_groupthreshold = 0.
    32         self.deviatoricerror_maximum = 0.
    33         self.restart = 0.
    34     #set defaults
     17        self.err = 0
     18        self.keepmetric = 0
     19        self.gradation = 0
     20        self.groundingline_resolution = 0
     21        self.groundingline_distance = 0
     22        self.icefront_resolution = 0
     23        self.icefront_distance = 0
     24        self.thicknesserror_resolution = 0
     25        self.thicknesserror_threshold = 0
     26        self.thicknesserror_groupthreshold = 0
     27        self.thicknesserror_maximum = 0
     28        self.deviatoricerror_resolution = 0
     29        self.deviatoricerror_threshold = 0
     30        self.deviatoricerror_groupthreshold = 0
     31        self.deviatoricerror_maximum = 0
     32        self.restart = 0
     33
    3534        self.setdefaultparameters()
    3635    #}}}
    3736
    38     def __repr__(self):  # {{{
    39         string = "   amr parameters:"
    40         string = "%s\n%s" % (string, fielddisplay(self, "hmin", "minimum element length"))
    41         string = "%s\n%s" % (string, fielddisplay(self, "hmax", "maximum element length"))
    42         string = "%s\n%s" % (string, fielddisplay(self, "fieldname", "name of input that will be used to compute the metric (should be an input of FemModel)"))
    43         string = "%s\n%s" % (string, fielddisplay(self, "keepmetric", "indicates whether the metric should be kept every remeshing time"))
    44         string = "%s\n%s" % (string, fielddisplay(self, "gradation", "maximum ratio between two adjacent edges"))
    45         string = "%s\n%s" % (string, fielddisplay(self, "groundingline_resolution", "element length near the grounding line"))
    46         string = "%s\n%s" % (string, fielddisplay(self, "groundingline_distance", "distance around the grounding line which elements will be refined"))
    47         string = "%s\n%s" % (string, fielddisplay(self, "icefront_resolution", "element length near the ice front"))
    48         string = "%s\n%s" % (string, fielddisplay(self, "icefront_distance", "distance around the ice front which elements will be refined"))
    49         string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_resolution", "element length when thickness error estimator is used"))
    50         string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_threshold", "maximum threshold thickness error permitted"))
    51         string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_groupthreshold", "maximum group threshold thickness error permitted"))
    52         string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_maximum", "maximum thickness error permitted"))
    53         string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_resolution", "element length when deviatoric stress error estimator is used"))
    54         string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_threshold", "maximum threshold deviatoricstress error permitted"))
    55         string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_groupthreshold", "maximum group threshold deviatoric stress error permitted"))
    56         string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_maximum", "maximum deviatoricstress error permitted"))
    57         string = "%s\n%s" % (string, fielddisplay(self, "restart", "indicates if ReMesh() will call before first time step"))
    58         return string
     37    def __repr__(self): #{{{
     38        s = '   amr parameters:\n'
     39        s += '{}\n'.format(fielddisplay(self, 'hmin', 'minimum element length'))
     40        s += '{}\n'.format(fielddisplay(self, 'hmax', 'maximum element length'))
     41        s += '{}\n'.format(fielddisplay(self, 'fieldname', 'name of input that will be used to compute the metric (should be an input of FemModel)'))
     42        s += '{}\n'.format(fielddisplay(self, 'keepmetric', 'indicates whether the metric should be kept every remeshing time'))
     43        s += '{}\n'.format(fielddisplay(self, 'gradation', 'maximum ratio between two adjacent edges'))
     44        s += '{}\n'.format(fielddisplay(self, 'groundingline_resolution', 'element length near the grounding line'))
     45        s += '{}\n'.format(fielddisplay(self, 'groundingline_distance', 'distance around the grounding line which elements will be refined'))
     46        s += '{}\n'.format(fielddisplay(self, 'icefront_resolution', 'element length near the ice front'))
     47        s += '{}\n'.format(fielddisplay(self, 'icefront_distance', 'distance around the ice front which elements will be refined'))
     48        s += '{}\n'.format(fielddisplay(self, 'thicknesserror_resolution', 'element length when thickness error estimator is used'))
     49        s += '{}\n'.format(fielddisplay(self, 'thicknesserror_threshold', 'maximum threshold thickness error permitted'))
     50        s += '{}\n'.format(fielddisplay(self, 'thicknesserror_groupthreshold', 'maximum group threshold thickness error permitted'))
     51        s += '{}\n'.format(fielddisplay(self, 'thicknesserror_maximum', 'maximum thickness error permitted'))
     52        s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_resolution', 'element length when deviatoric stress error estimator is used'))
     53        s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_threshold', 'maximum threshold deviatoricstress error permitted'))
     54        s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_groupthreshold', 'maximum group threshold deviatoric stress error permitted'))
     55        s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_maximum', 'maximum deviatoricstress error permitted'))
     56        s += '{}\n'.format(fielddisplay(self, 'restart', 'indicates if ReMesh() will call before first time step'))
     57        return s
    5958    #}}}
    6059
    61     def setdefaultparameters(self):  # {{{
    62         self.hmin = 100.
    63         self.hmax = 100.e3
     60    def setdefaultparameters(self): #{{{
     61        self.hmin = 100
     62        self.hmax = 100e3
     63
     64        # Fields
    6465        self.fieldname = 'Vel'
    65         self.err = 3.
     66        self.err = 3
     67
     68        # Keep metric?
    6669        self.keepmetric = 1
     70
     71        # Control of element lengths
    6772        self.gradation = 1.5
    68         self.groundingline_resolution = 500.
     73
     74        # Other criteria
     75        self.groundingline_resolution = 500
    6976        self.groundingline_distance = 0
    70         self.icefront_resolution = 500.
     77        self.icefront_resolution = 500
    7178        self.icefront_distance = 0
    72         self.thicknesserror_resolution = 500.
     79        self.thicknesserror_resolution = 500
    7380        self.thicknesserror_threshold = 0
    7481        self.thicknesserror_groupthreshold = 0
    7582        self.thicknesserror_maximum = 0
    76         self.deviatoricerror_resolution = 500.
     83        self.deviatoricerror_resolution = 500
    7784        self.deviatoricerror_threshold = 0
    7885        self.deviatoricerror_groupthreshold = 0
    7986        self.deviatoricerror_maximum = 0
    80         self.restart = 0.
     87
     88        # Is restart? This calls femmodel->ReMesh() before first time step.
     89        self.restart = 0
    8190        return self
    8291    #}}}
    8392
    84     def checkconsistency(self, md, solution, analyses):  # {{{
     93    def checkconsistency(self, md, solution, analyses): #{{{
    8594        md = checkfield(md, 'fieldname', 'amr.hmax', 'numel', [1], '>', 0, 'NaN', 1)
    8695        md = checkfield(md, 'fieldname', 'amr.hmin', 'numel', [1], '>', 0, '<', self.hmax, 'NaN', 1)
    … …  
    101110        md = checkfield(md, 'fieldname', 'amr.restart', 'numel', [1], '>=', 0, '<=', 1, 'NaN', 1)
    102111        return md
    103     # }}}
     112   # }}}
    104113
    105     def marshall(self, prefix, md, fid):  # {{{
     114    def marshall(self, prefix, md, fid): #{{{
    106115        WriteData(fid, prefix, 'name', 'md.amr.type', 'data', 1, 'format', 'Integer')
    107116        WriteData(fid, prefix, 'object', self, 'fieldname', 'hmin', 'format', 'Double')
    … …  
    124133        WriteData(fid, prefix, 'object', self, 'fieldname', 'deviatoricerror_maximum', 'format', 'Double')
    125134        WriteData(fid, prefix, 'object', self, 'class', 'amr', 'fieldname', 'restart', 'format', 'Integer')
    126     # }}}
     135    #}}}

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

    r26231 r26301  
    1313        end
    1414        methods
    15                 function self = extrude(self,md) % {{{
    16                         self.sea_surface_height_above_geoid=project3d(md,'vector',self.sea_surface_height_above_geoid,'type','node','layer',1);
    17                         self.sea_water_pressure_at_sea_floor=project3d(md,'vector',self.sea_water_pressure_at_sea_floor,'type','node','layer',1);
    18                 end % }}}
    1915                function self = dsl(varargin) % {{{
    2016                        switch nargin
    … …  
    3430
    3531                end % }}}
     32                function disp(self) % {{{
     33
     34                        disp(sprintf('   dsl parameters:'));
     35                        fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).');
     36                        fielddisplay(self,'sea_surface_height_above_geoid','Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).');
     37                        fielddisplay(self,'sea_water_pressure_at_sea_floor','Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).');
     38
     39                end % }}}
    3640                function md = checkconsistency(self,md,solution,analyses) % {{{
    3741
    … …  
    4953                       
    5054                end % }}}
    51                 function disp(self) % {{{
    52 
    53                         disp(sprintf('   dsl parameters:'));
    54                         fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).');
    55                         fielddisplay(self,'sea_surface_height_above_geoid','Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).');
    56                         fielddisplay(self,'sea_water_pressure_at_sea_floor','Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).');
     55                function marshall(self,prefix,md,fid) % {{{
     56                        yts=md.constants.yts;
     57                        WriteData(fid,prefix,'name','md.dsl.model','data',1,'format','Integer');
     58                        WriteData(fid,prefix,'object',self,'fieldname','global_average_thermosteric_sea_level','format','DoubleMat','mattype',2,'timeseries',1,'timeserieslength',2,'yts',yts); %mattype 2, because we are sending a GMSL value identical everywhere on each element.
     59                        WriteData(fid,prefix,'object',self,'fieldname','sea_surface_height_above_geoid','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); %mattype 1 because we specify DSL at vertex locations.
     60                        WriteData(fid,prefix,'object',self,'fieldname','sea_water_pressure_at_sea_floor','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); %mattype 1 because we specify bottom pressure at vertex locations.
    5761
    5862                end % }}}
    59                 function marshall(self,prefix,md,fid) % {{{
    60 
    61                         WriteData(fid,prefix,'name','md.dsl.model','data',1,'format','Integer');
    62                         WriteData(fid,prefix,'object',self,'fieldname','global_average_thermosteric_sea_level','format','DoubleMat','mattype',2,'timeseries',1,'timeserieslength',2,'yts',md.constants.yts); %mattype 2, because we are sending a GMSL value identical everywhere on each element.
    63                         WriteData(fid,prefix,'object',self,'fieldname','sea_surface_height_above_geoid','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); %mattype 1 because we specify DSL at vertex locations.
    64                         WriteData(fid,prefix,'object',self,'fieldname','sea_water_pressure_at_sea_floor','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); %mattype 1 because we specify bottom pressure at vertex locations.
    65 
    66                 end % }}}
    67                 function savemodeljs(self,fid,modelname) % {{{
    68 
    69                         writejs1Darray(fid,[modelname '.dsl.global_average_thermosteric_sea_level'],self.global_average_thermosteric_sea_level);
    70                         writejs1Darray(fid,[modelname '.dsl.sea_surface_height_above_geoid'],self.sea_surface_height_above_geoid);
    71                         writejs1Darray(fid,[modelname '.dsl.sea_water_pressure_at_sea_floor'],self.sea_water_pressure_at_sea_floor);
    72 
     63                function self = extrude(self,md) % {{{
     64                        self.sea_surface_height_above_geoid=project3d(md,'vector',self.sea_surface_height_above_geoid,'type','node','layer',1);
     65                        self.sea_water_pressure_at_sea_floor=project3d(md,'vector',self.sea_water_pressure_at_sea_floor,'type','node','layer',1);
    7366                end % }}}
    7467                function self = initialize(self,md) % {{{
    … …  
    8780                        end
    8881                end % }}}
    89        
     82                function savemodeljs(self,fid,modelname) % {{{
     83
     84                        writejs1Darray(fid,[modelname '.dsl.global_average_thermosteric_sea_level'],self.global_average_thermosteric_sea_level);
     85                        writejs1Darray(fid,[modelname '.dsl.sea_surface_height_above_geoid'],self.sea_surface_height_above_geoid);
     86                        writejs1Darray(fid,[modelname '.dsl.sea_water_pressure_at_sea_floor'],self.sea_water_pressure_at_sea_floor);
     87
     88                end % }}}
    9089        end
    9190end

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

    r26066 r26301  
    88
    99class dsl(object):
    10     """DSL - slass definition
     10    """DSL - class definition
    1111
    1212    Usage:
    13         dsl = dsl()
     13        dsl = dsl() # dynamic sea level class, based on CMIP5 outputs
    1414    """
    1515
    16     def __init__(self): #{{{
     16    def __init__(self, *args): #{{{
    1717        self.global_average_thermosteric_sea_level = np.nan # Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).
    1818        self.sea_surface_height_above_geoid        = np.nan # Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).
    1919        self.sea_water_pressure_at_sea_floor       = np.nan # Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).
     20
     21        if len(args) == 0:
     22            self.setdefaultparameters()
     23        else:
     24            raise Exception('constructor not supported')
    2025    #}}}
    2126
    … …  
    2833    #}}}
    2934
    30     def defaultoutputs(self, md): #{{{
    31         return []
    32     #}}}
    33 
    34     def initialize(self, md): #{{{
    35         if np.isnan(self.global_average_thermosteric_sea_level):
    36             self.global_average_thermosteric_sea_level = np.array([0, 0]).reshape(-1, 1)
    37             print('      no dsl.global_average_thermosteric_sea_level specified: transient values set at zero')
    38 
    39         if np.isnan(self.sea_surface_height_above_geoid):
    40             self.sea_surface_height_above_geoid = np.array(np.zeros(md.mesh.numberofvertices)).reshape(-1, 1)
    41             print('      no dsl.sea_surface_height_above_geoid specified: transient values set at zero')
    42 
    43         if np.isnan(self.sea_water_pressure_at_sea_floor):
    44             self.sea_water_pressure_at_sea_floor = np.array(np.zeros(md.mesh.numberofvertices)).reshape(-1, 1)
    45             print('      no dsl.sea_water_pressure_at_sea_floor specified: transient values set at zero')
     35    def setdefaultparameters(self): #{{{
     36        self.global_average_thermosteric_sea_level = np.nan
     37        self.sea_surface_height_above_geoid        = np.nan
     38        self.sea_water_pressure_at_sea_floor       = np.nan
    4639    #}}}
    4740
    4841    def checkconsistency(self, md, solution, analyses): #{{{
    4942        # Early return
    50         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
     43        if ('SealevelchangeAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
    5144            return md
    5245
    … …  
    6154    # }}}
    6255
     56    def marshall(self, prefix, md, fid):   #{{{
     57        yts = md.constants.yts
     58        WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer')
     59        WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'DoubleMat', 'mattype', 2, 'timeseries', 1, 'yts', yts) # mattype 2, because we are sending a GMSL value identical everywhere on each element.
     60        WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) # mattype 1 because we specify DSL at vertex locations.
     61        WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) # mattype 1 because we specify bottom pressure at vertex locations.
     62    # }}}
     63
    6364    def extrude(self, md): #{{{
    64         self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node')
    65         self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node')
     65        self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node', 'layer', 1)
     66        self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node', 'layer', 1)
    6667        return self
    6768    #}}}
    6869
    69     def marshall(self, prefix, md, fid):   #{{{
    70         yts = md.constants.yts
    71         WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer')
    72         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', 1+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) # mattype 2, because we are sending a GMSL value identical everywhere on each element.
    73         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) # mattype 1 because we specify DSL at vertex locations.
    74         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_surface_height_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts) # mattype 1 because we specify bottom pressure at vertex locations.
    75     # }}}
     70
     71    def initialize(self, md): #{{{
     72        if np.isnan(self.global_average_thermosteric_sea_level):
     73            self.global_average_thermosteric_sea_level = np.array([0, 0]).reshape(-1, 1)
     74            print('      no dsl.global_average_thermosteric_sea_level specified: transient values set at zero')
     75
     76        if np.isnan(self.sea_surface_height_above_geoid):
     77            self.sea_surface_height_above_geoid = np.append(np.zeros((md.mesh.numberofvertices, 1)), 0).reshape(-1, 1)
     78            print('      no dsl.sea_surface_height_above_geoid specified: transient values set at zero')
     79
     80        if np.isnan(self.sea_water_pressure_at_sea_floor):
     81            self.sea_water_pressure_at_sea_floor = np.append(np.zeros((md.mesh.numberofvertices, 1)), 0).reshape(-1, 1)
     82            print('      no dsl.sea_water_pressure_at_sea_floor specified: transient values set at zero')
     83    #}}}

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

    r26242 r26301  
    66classdef fourierlove
    77        properties (SetAccess=public)
    8                 nfreq                      = 0;
    9                 frequencies                = 0;
    10                 sh_nmax                    = 0;
    11                 sh_nmin                    = 0;
    12                 g0                         = 0;
    13                 r0                         = 0;
    14                 mu0                        = 0;
    15                 Gravitational_Constant     = 0;
    16                 allow_layer_deletion       = 0;
    17                 underflow_tol              = 0;
    18                 integration_steps_per_layer= 0;
    19                 istemporal                 = 0;
    20                 n_temporal_iterations      = 0;
    21                 time                               = 0;
    22                 love_kernels               = 0;
    23                 forcing_type               = 0;
    24                 inner_core_boundary            = 0;
    25                 core_mantle_boundary       = 0;
    26                 complex_computation        = 0;
     8                nfreq                       = 0;
     9                frequencies                 = 0;
     10                sh_nmax                     = 0;
     11                sh_nmin                     = 0;
     12                g0                          = 0;
     13                r0                          = 0;
     14                mu0                         = 0;
     15                Gravitational_Constant      = 0;
     16                allow_layer_deletion        = 0;
     17                underflow_tol               = 0;
     18                integration_steps_per_layer = 0;
     19                istemporal                  = 0;
     20                n_temporal_iterations       = 0;
     21                time                        = 0;
     22                love_kernels                = 0;
     23                forcing_type                = 0;
     24                inner_core_boundary         = 0;
     25                core_mantle_boundary        = 0;
     26                complex_computation         = 0;
    2727
    2828        end
    … …  
    3434        end
    3535        methods
    36                 function self = extrude(self,md) % {{{
    37                 end % }}}
    3836                function self = fourierlove(varargin) % {{{
    3937                        switch nargin
    … …  
    5351                        self.g0=9.81; % m/s^2;
    5452                        self.r0=6371*1e3; %m;
    55                         self.mu0=10^11; % Pa
     53                        self.mu0=1e11; % Pa
    5654                        self.Gravitational_Constant=6.67259e-11; % m^3 kg^-1 s^-2
    5755                        self.allow_layer_deletion=1;
    … …  
    6866                end % }}}
    6967                function disp(self) % {{{
    70                         fielddisplay(self,'nfreq','number of frequencies sampled (default 1, elastic) [Hz]');
     68                        fielddisplay(self,'nfreq','number of frequencies sampled (default: 1, elastic) [Hz]');
    7169                        fielddisplay(self,'frequencies','frequencies sampled (convention defaults to 0 for the elastic case) [Hz]');
    72                         fielddisplay(self,'sh_nmax','maximum spherical harmonic degree (default 256, .35 deg, or 40 km at equator)');
    73                         fielddisplay(self,'sh_nmin','minimum spherical harmonic degree (default 1)');
    74                         fielddisplay(self,'g0','adimensioning constant for gravity (default 10) [m/s^2]');
    75                         fielddisplay(self,'r0','adimensioning constant for radius (default 6371*10^3) [m]');
    76                         fielddisplay(self,'mu0','adimensioning constant for stress (default 10^11) [Pa]');
    77                         fielddisplay(self,'Gravitational_Constant','Newtonian constant of gravitation (default 6.67259e-11 [m^3 kg^-1 s^-2])');
    78                         fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default 1)');
    79                         fielddisplay(self,'underflow_tol','threshold of deep to surface love number ratio to trigger the deletion of layers (default 2.2204460492503131E-016)');
    80                         fielddisplay(self,'integration_steps_per_layer','number of radial steps to propagate the yi system from the bottom to the top of each layer (default 100)');
    81                         fielddisplay(self,'istemporal',{'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'});
    82                         fielddisplay(self,'n_temporal_iterations','max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default 8)');
    83                         fielddisplay(self,'time','time vector for deformation if istemporal (default 0) [s]');
    84                         fielddisplay(self,'love_kernels','compute love numbers at depth? (default 0)');
    85                         fielddisplay(self,'forcing_type',{'integer indicating the nature and depth of the forcing for the Love number calculation (default 11) :','1:  Inner core boundary -- Volumic Potential','2:  Inner core boundary -- Pressure','3:  Inner core boundary -- Loading','4:  Inner core boundary -- Tangential traction','5:  Core mantle boundary -- Volumic Potential','6:  Core mantle boundary -- Pressure','7:  Core mantle boundary -- Loading','8:  Core mantle boundary -- Tangential traction','9:  Surface -- Volumic Potential','10: Surface -- Pressure','11: Surface -- Loading','12: Surface -- Tangential traction '});
    86                         fielddisplay(self,'inner_core_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default 1)');
    87                         fielddisplay(self,'core_mantle_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default 2)');
     70                        fielddisplay(self,'sh_nmax','maximum spherical harmonic degree (default: 256, .35 deg, or 40 km at equator)');
     71                        fielddisplay(self,'sh_nmin','minimum spherical harmonic degree (default: 1)');
     72                        fielddisplay(self,'g0','adimensioning constant for gravity (default: 10) [m/s^2]');
     73                        fielddisplay(self,'r0','adimensioning constant for radius (default: 6371*10^3) [m]');
     74                        fielddisplay(self,'mu0','adimensioning constant for stress (default: 10^11) [Pa]');
     75                        fielddisplay(self,'Gravitational_Constant','Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])');
     76                        fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)');
     77                        fielddisplay(self,'underflow_tol','threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)');
     78                        fielddisplay(self,'integration_steps_per_layer','number of radial steps to propagate the yi system from the bottom to the top of each layer (default: 100)');
     79                        fielddisplay(self,'istemporal',{'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default: 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'});
     80                        fielddisplay(self,'n_temporal_iterations','max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default: 8)');
     81                        fielddisplay(self,'time','time vector for deformation if istemporal (default: 0) [s]');
     82                        fielddisplay(self,'love_kernels','compute love numbers at depth? (default: 0)');
     83                        fielddisplay(self,'forcing_type',{'integer indicating the nature and depth of the forcing for the Love number calculation (default: 11):','1:  Inner core boundary -- Volumic Potential','2:  Inner core boundary -- Pressure','3:  Inner core boundary -- Loading','4:  Inner core boundary -- Tangential traction','5:  Core mantle boundary -- Volumic Potential','6:  Core mantle boundary -- Pressure','7:  Core mantle boundary -- Loading','8:  Core mantle boundary -- Tangential traction','9:  Surface -- Volumic Potential','10: Surface -- Pressure','11: Surface -- Loading','12: Surface -- Tangential traction '});
     84                        fielddisplay(self,'inner_core_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default: 1)');
     85                        fielddisplay(self,'core_mantle_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default: 2)');
    8886
    8987                end % }}}
    … …  
    112110                                md = checkfield(md,'fieldname','love.time','NaN',1,'Inf',1,'numel',md.love.nfreq/2/md.love.n_temporal_iterations);
    113111                        end
    114                         if md.love.sh_nmin<=1 & (md.love.forcing_type==9 || md.love.forcing_type==5 || md.love.forcing_type==1)
    115                                 error('Degree 1 not supported for Volumetric Potential forcing. Use sh_min>=2 for this kind of calculation.')
     112                        if md.love.sh_nmin<=1 & (md.love.forcing_type==1 || md.love.forcing_type==5 || md.love.forcing_type==9)
     113                                error(['Degree 1 not supported for forcing type ' num2str(md.love.forcing_type) '. Use sh_min>=2 for this kind of calculation.'])
    116114                        end
    117115
    … …  
    130128                end % }}}
    131129                function marshall(self,prefix,md,fid) % {{{
    132                
     130
    133131                        WriteData(fid,prefix,'object',self,'fieldname','nfreq','format','Integer');
    134132                        WriteData(fid,prefix,'object',self,'fieldname','frequencies','format','DoubleMat','mattype',3);
    … …  
    152150
    153151                end % }}}
     152                function self = extrude(self,md) % {{{
     153                end % }}}
    154154                function savemodeljs(self,fid,modelname) % {{{
    155155                        error('not implemented yet!');

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

    r26059 r26301  
     1import numpy as np
     2
    13from fielddisplay import fielddisplay
    24from checkfield import checkfield
    … …  
    57
    68class fourierlove(object):
    7     """FOURIERLOVE - Fourier Love Number class definition
     9    """FOURIERLOVE - class definition
    810
    911    Usage:
    10         fourierlove = fourierlove()
     12        md.love = fourierlove()
    1113    """
    12     def __init__(self):  # {{{
     14
     15    def __init__(self): #{{{
    1316        self.nfreq = 0
    1417        self.frequencies = 0
    … …  
    1821        self.r0 = 0
    1922        self.mu0 = 0
     23        self.Gravitational_Constant = 0
    2024        self.allow_layer_deletion = 0
     25        self.underflow_tol = 0
     26        self.integration_steps_per_layer = 0
     27        self.istemporal = 0
     28        self.n_temporal_iterations = 0
     29        self.time = 0
    2130        self.love_kernels = 0
    2231        self.forcing_type = 0
     32        self.inner_core_boundary = 0
     33        self.core_mantle_boundary = 0
     34        self.complex_computation = 0
    2335
    2436        self.setdefaultparameters()
    2537    #}}}
    2638
    27     def __repr__(self):  # {{{
    28         # TODO:
    29         # - Convert all formatting to calls to <string>.format (see any
    30         #   already converted <class>.__repr__ method for examples)
    31         #
    32         string = '   Fourier Love class:'
    33         string = "%s\n%s" % (string, fielddisplay(self, 'nfreq', 'number of frequencies sampled (default 1, elastic) [Hz]'))
    34         string = "%s\n%s" % (string, fielddisplay(self, 'frequencies', 'frequencies sampled (convention defaults to 0 for the elastic case) [Hz]'))
    35         string = "%s\n%s" % (string, fielddisplay(self, 'sh_nmax', 'maximum spherical harmonic degree (default 256, .35 deg, or 40 km at equator)'))
    36         string = "%s\n%s" % (string, fielddisplay(self, 'sh_nmin', 'minimum spherical harmonic degree (default 1)'))
    37         string = "%s\n%s" % (string, fielddisplay(self, 'g0', 'adimensioning constant for gravity (default 10) [m / s^2]'))
    38         string = "%s\n%s" % (string, fielddisplay(self, 'r0', 'adimensioning constant for radius (default 6378e3) [m]'))
    39         string = "%s\n%s" % (string, fielddisplay(self, 'mu0', 'adimensioning constant for stress (default 1.0e11) [Pa]'))
    40         string = "%s\n%s" % (string, fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default 1)'))
    41         string = "%s\n%s" % (string, fielddisplay(self, 'love_kernels', 'compute love numbers at depth? (default 0)'))
    42         string = "%s\n%s" % (string, fielddisplay(self, 'forcing_type', 'integer indicating the nature and depth of the forcing for the Love number calculation (default 11) :'))
    43         string = "%s\n%s" % (string, '                                                     1:  Inner core boundary -- Volumic Potential')
    44         string = "%s\n%s" % (string, '                                                     2:  Inner core boundary --  Pressure')
    45         string = "%s\n%s" % (string, '                                                     3:  Inner core boundary --  Loading')
    46         string = "%s\n%s" % (string, '                                                     4:  Inner core boundary --  Tangential traction')
    47         string = "%s\n%s" % (string, '                                                     5:  Core mantle boundary --  Volumic Potential')
    48         string = "%s\n%s" % (string, '                                                     6:  Core mantle boundary --  Pressure')
    49         string = "%s\n%s" % (string, '                                                     7:  Core mantle boundary --  Loading')
    50         string = "%s\n%s" % (string, '                                                     8:  Core mantle boundary --  Tangential traction')
    51         string = "%s\n%s" % (string, '                                                     9:  Surface--  Volumic Potential')
    52         string = "%s\n%s" % (string, '                                                     10: Surface--  Pressure')
    53         string = "%s\n%s" % (string, '                                                     11: Surface--  Loading')
    54         string = "%s\n%s" % (string, '                                                     12: Surface--  Tangential traction ')
     39    def __repr__(self): #{{{
     40        s = '   Fourier Love class:\n'
     41        s += '{}\n'.format(fielddisplay(self, 'nfreq', 'number of frequencies sampled (default: 1, elastic) [Hz]'))
     42        s += '{}\n'.format(fielddisplay(self, 'frequencies', 'frequencies sampled (convention defaults to 0 for the elastic case) [Hz]'))
     43        s += '{}\n'.format(fielddisplay(self, 'sh_nmax', 'maximum spherical harmonic degree (default: 256, .35 deg, or 40 km at equator)'))
     44        s += '{}\n'.format(fielddisplay(self, 'sh_nmin', 'minimum spherical harmonic degree (default: 1)'))
     45        s += '{}\n'.format(fielddisplay(self, 'g0', 'adimensioning constant for gravity (default: 10) [m / s^2]'))
     46        s += '{}\n'.format(fielddisplay(self, 'r0', 'adimensioning constant for radius (default: 6378e3) [m]'))
     47        s += '{}\n'.format(fielddisplay(self, 'mu0', 'adimensioning constant for stress (default: 1.0e11) [Pa]'))
     48        s += '{}\n'.format(fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)'))
     49        s += '{}\n'.format(fielddisplay(self, 'Gravitational_Constant', 'Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])'))
     50        s += '{}\n'.format(fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)'))
     51        s += '{}\n'.format(fielddisplay(self, 'underflow_tol', 'threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)'))
     52        s += '{}\n'.format(fielddisplay(self, 'integration_steps_per_layer', 'number of radial steps to propagate the yi system from the bottom to the top of each layer (default: 100)'))
     53        s += '{}\n'.format(fielddisplay(self, 'istemporal', {'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default: 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'}))
     54        s += '{}\n'.format(fielddisplay(self, 'n_temporal_iterations', 'max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default: 8)'))
     55        s += '{}\n'.format(fielddisplay(self, 'time', 'time vector for deformation if istemporal (default: 0) [s]'))
     56        s += '{}\n'.format(fielddisplay(self, 'love_kernels', 'compute love numbers at depth? (default: 0)'))
     57        s += '{}\n'.format(fielddisplay(self, 'forcing_type', 'integer indicating the nature and depth of the forcing for the Love number calculation (default: 11):'))
     58        s += '{}\n'.format('                                                     1:  Inner core boundary -- Volumic Potential')
     59        s += '{}\n'.format('                                                     2:  Inner core boundary --  Pressure')
     60        s += '{}\n'.format('                                                     3:  Inner core boundary --  Loading')
     61        s += '{}\n'.format('                                                     4:  Inner core boundary --  Tangential traction')
     62        s += '{}\n'.format('                                                     5:  Core mantle boundary --  Volumic Potential')
     63        s += '{}\n'.format('                                                     6:  Core mantle boundary --  Pressure')
     64        s += '{}\n'.format('                                                     7:  Core mantle boundary --  Loading')
     65        s += '{}\n'.format('                                                     8:  Core mantle boundary --  Tangential traction')
     66        s += '{}\n'.format('                                                     9:  Surface--  Volumic Potential')
     67        s += '{}\n'.format('                                                     10: Surface--  Pressure')
     68        s += '{}\n'.format('                                                     11: Surface--  Loading')
     69        s += '{}\n'.format('                                                     12: Surface--  Tangential traction ')
     70        s += '{}\n'.format(fielddisplay(self, 'inner_core_boundary', 'interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default: 1)'))
     71        s += '{}\n'.format(fielddisplay(self, 'core_mantle_boundary', 'interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default: 2)'))
    5572
    56         return string
     73        return s
    5774    #}}}
    5875
    59     def extrude(self, md):  # {{{
    60         return self
     76    def setdefaultparameters(self): #{{{
     77        # We setup an elastic love number computation by default
     78        self.nfreq = 1
     79        self.frequencies = [0] # Hz
     80        self.sh_nmax = 256 # .35 degree, 40 km at the equator
     81        self.sh_nmin = 1
     82        # Work on matlab script for computing g0 for given Earth's structure
     83        self.g0 = 9.81 # m/s^2
     84        self.r0 = 6371 * 1e3 # m
     85        self.mu0 = 1e11 # Pa
     86        self.Gravitational_Constant = 6.67259e-11 # m^3 kg^-1 s^-2
     87        self.allow_layer_deletion = 1
     88        self.underflow_tol = 1e-16 # Threshold of deep to surface love number ratio to trigger the deletion of layer
     89        self.integration_steps_per_layer = 100
     90        self.istemporal = 0
     91        self.n_temporal_iterations = 8
     92        self.time = [0] # s
     93        self.love_kernels = 0
     94        self.forcing_type = 11 # Surface loading
     95        self.inner_core_boundary = 1
     96        self.core_mantle_boundary = 2
     97        self.complex_computation = 0
    6198    #}}}
    6299
    63     def setdefaultparameters(self):  # {{{
    64         #we setup an elastic love number computation by default.
    65         self.nfreq = 1
    66         self.frequencies = [0]  #Hz
    67         self.sh_nmax = 256  # .35 degree, 40 km at the equator.
    68         self.sh_nmin = 1
    69         self.g0 = 9.81 # m/s^2
    70         self.r0 = 6371e3 #m
    71         self.mu0 = 1e11 # Pa
    72         self.allow_layer_deletion = 1
    73         self.love_kernels = 0
    74         self.forcing_type = 11
     100    def checkconsistency(self, md, solution, analyses): #{{{
     101        if 'LoveAnalysis' not in analyses:
     102            return md
    75103
    76         return self
     104        md = checkfield(md, 'fieldname', 'love.nfreq', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     105        md = checkfield(md, 'fieldname', 'love.frequencies', 'NaN', 1, 'Inf', 1, 'numel', md.love.nfreq)
     106        md = checkfield(md, 'fieldname', 'love.sh_nmax', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     107        md = checkfield(md, 'fieldname', 'love.sh_nmin', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     108        md = checkfield(md, 'fieldname', 'love.g0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     109        md = checkfield(md, 'fieldname', 'love.r0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     110        md = checkfield(md, 'fieldname', 'love.mu0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     111        md = checkfield(md, 'fieldname', 'love.Gravitational_Constant', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     112        md = checkfield(md, 'fieldname', 'love.allow_layer_deletion', 'values', [0, 1])
     113        md = checkfield(md, 'fieldname', 'love.underflow_tol', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     114        md = checkfield(md, 'fieldname', 'love.integration_steps_per_layer', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     115        md = checkfield(md, 'fieldname', 'love.love_kernels', 'values', [0, 1])
     116        md = checkfield(md, 'fieldname', 'love.forcing_type', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', 12)
     117        md = checkfield(md, 'fieldname', 'love.complex_computation', 'NaN', 1, 'Inf', 1, 'numel', 1, 'values', [0, 1])
     118
     119        md = checkfield(md, 'fieldname', 'love.istemporal', 'values', [0, 1])
     120        if md.love.istemporal:
     121            md = checkfield(md, 'fieldname', 'love.n_temporal_iterations', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     122            md = checkfield(md, 'fieldname', 'love.time', 'NaN', 1, 'Inf', 1, 'numel', md.love.nfreq / 2 / md.love.n_temporal_iterations)
     123        if md.love.sh_nmin <= 1 and (md.love.forcing_type == 1 or md.love.forcing_type == 5 or md.love.forcing_type == 9):
     124            raise RuntimeError('Degree 1 not supported for forcing type {}. Use sh_min >= 2 for this kind of calculation.'.format(md.love.forcing_type))
     125
     126        # Need 'litho' material
     127        if not hasattr(md.materials, 'materials') or 'litho' not in md.materials.nature:
     128            raise RuntimeError('Need a \'litho\' material to run a Fourier Love number analysis')
     129
     130        mat = np.where(np.array(nature) == 'litho')[0]
     131        if md.love.forcing_type <= 4:
     132            md = checkfield(md, 'fieldname', 'love.inner_core_boundary', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', md.materials[mat].numlayers)
     133        elif md.love.forcing_type <= 8:
     134            md = checkfield(md, 'fieldname', 'love.core_mantle_boundary', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', md.materials[mat].numlayers)
     135
     136        return md
    77137    #}}}
    78138
    79     def checkconsistency(self, md, solution, analyses):  # {{{
    80         if 'LoveAnalysis' not in analyses:
    81             return md
    82         md = checkfield(md, 'fieldname', 'love.nfreq', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    83         md = checkfield(md, 'fieldname', 'love.frequencies', 'NaN', 1, 'Inf', 1, 'numel', [md.love.nfreq])
    84         md = checkfield(md, 'fieldname', 'love.sh_nmax', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    85         md = checkfield(md, 'fieldname', 'love.sh_nmin', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    86         md = checkfield(md, 'fieldname', 'love.g0', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    87         md = checkfield(md, 'fieldname', 'love.r0', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    88         md = checkfield(md, 'fieldname', 'love.mu0', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    89         md = checkfield(md, 'fieldname', 'love.allow_layer_deletion', 'values', [0, 1])
    90         md = checkfield(md, 'fieldname', 'love.love_kernels', 'values', [0, 1])
    91         md = checkfield(md, 'fieldname', 'love.forcing_type', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0, '<=', 12)
    92         if md.love.sh_nmin <= 1 and md.love.forcing_type == 9:
    93             raise RuntimeError("Degree 1 not supported for Volumetric Potential forcing. Use sh_min >= 2 for this kind of calculation.")
    94 
    95         # need 'litho' material
    96         if not hasattr(md.materials, 'materials') or 'litho' not in md.materials.nature:
    97             raise RuntimeError('Need a \'litho\' material to run a Fourier Love number analysis')
    98         return md
    99     # }}}
    100 
    101     def marshall(self, prefix, md, fid):  # {{{
     139    def marshall(self, prefix, md, fid): #{{{
    102140        WriteData(fid, prefix, 'object', self, 'fieldname', 'nfreq', 'format', 'Integer')
    103         WriteData(fid, prefix, 'object', self, 'fieldname', 'frequencies', 'format', 'DoubleMat', 'mattype', 3)
     141        WriteData(fid, prefix, 'object', self, 'fieldname', 'frequencies', 'format', 'DoubleMat', 'mattype',3)
    104142        WriteData(fid, prefix, 'object', self, 'fieldname', 'sh_nmax', 'format', 'Integer')
    105143        WriteData(fid, prefix, 'object', self, 'fieldname', 'sh_nmin', 'format', 'Integer')
    … …  
    107145        WriteData(fid, prefix, 'object', self, 'fieldname', 'r0', 'format', 'Double')
    108146        WriteData(fid, prefix, 'object', self, 'fieldname', 'mu0', 'format', 'Double')
     147        WriteData(fid, prefix, 'object', self, 'fieldname', 'Gravitational_Constant', 'format', 'Double')
    109148        WriteData(fid, prefix, 'object', self, 'fieldname', 'allow_layer_deletion', 'format', 'Boolean')
     149        WriteData(fid, prefix, 'object', self, 'fieldname', 'underflow_tol', 'format', 'Double')
     150        WriteData(fid, prefix, 'object', self, 'fieldname', 'integration_steps_per_layer', 'format', 'Integer')
     151        WriteData(fid, prefix, 'object', self, 'fieldname', 'istemporal', 'format', 'Boolean')
     152        WriteData(fid, prefix, 'object', self, 'fieldname', 'n_temporal_iterations', 'format', 'Integer')
     153        WriteData(fid, prefix, 'object', self, 'fieldname', 'complex_computation', 'format', 'Boolean')
     154        # Note: no need to marshall the time vector, we have frequencies
    110155        WriteData(fid, prefix, 'object', self, 'fieldname', 'love_kernels', 'format', 'Boolean')
    111156        WriteData(fid, prefix, 'object', self, 'fieldname', 'forcing_type', 'format', 'Integer')
    112     # }}}
     157        WriteData(fid, prefix, 'object', self, 'fieldname', 'inner_core_boundary', 'format', 'Integer')
     158        WriteData(fid, prefix, 'object', self, 'fieldname', 'core_mantle_boundary', 'format', 'Integer')
     159    #}}}
     160
     161    def extrude(self, md): #{{{
     162        return self
     163    #}}}
     164
     165    def build_frequencies_from_time(self): #{{{
     166        if not self.istemporal:
     167            raise RuntimeError('cannot build frequencies for temporal love numbers if love.istemporal==0')
     168        print('Temporal love numbers: Overriding md.love.nfreq and md.love.frequencies')
     169        self.nfreq = len(self.time) * 2 * self.n_temporal_iterations
     170        self.frequencies = np.zeros((self.nfreq, 1))
     171        for i in range(len(self.time)):
     172            for j in range(2 * self.n_temporal_iterations):
     173                self.frequencies[(i - 1) * 2 * self.n_temporal_iterations + j] = j * np.log(2) / self.time[i] / 2 / np.pi
     174    #}}}

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

    r26299 r26301  
    3737
    3838    def setdefaultparameters(self): #{{{
    39         return self
     39        return
    4040    #}}}
    4141

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

    r26059 r26301  
    1111        end
    1212        methods
    13                 function self = extrude(self,md) % {{{
    14                 end % }}}
    1513                function self = hydrologyshreve(varargin) % {{{
    1614                        switch nargin
    … …  
    2321                        end
    2422                end % }}}
    25                 function list = defaultoutputs(self,md) % {{{
    26                         list = {'Watercolumn','HydrologyWaterVx','HydrologyWaterVy'};
    27                 end % }}}   
     23                function disp(self) % {{{
     24                        disp(sprintf('   hydrologyshreve solution parameters:'));
     25                        fielddisplay(self,'spcwatercolumn','water thickness constraints (NaN means no constraint) [m]');
     26                        fielddisplay(self,'stabilization','artificial diffusivity (default: 1). can be more than 1 to increase diffusivity.');
     27                        fielddisplay(self,'requested_outputs','additional outputs requested');
    2828
     29                end % }}}
    2930                function self = setdefaultparameters(self) % {{{
    3031
    … …  
    3435                end % }}}
    3536                function md = checkconsistency(self,md,solution,analyses) % {{{
    36                        
     37
    3738                        %Early return
    3839                        if ~ismember('HydrologyShreveAnalysis',analyses) |  (strcmp(solution,'TransientSolution') & md.transient.ishydrology==0), return; end
    … …  
    4142                        md = checkfield(md,'fieldname','hydrology.stabilization','>=',0);
    4243                end % }}}
    43                 function disp(self) % {{{
    44                         disp(sprintf('   hydrologyshreve solution parameters:'));
    45                         fielddisplay(self,'spcwatercolumn','water thickness constraints (NaN means no constraint) [m]');
    46                         fielddisplay(self,'stabilization','artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.');
    47       fielddisplay(self,'requested_outputs','additional outputs requested');
    48 
     44                function list = defaultoutputs(self,md) % {{{
     45                        list = {'Watercolumn','HydrologyWaterVx','HydrologyWaterVy'};
    4946                end % }}}
    5047                function marshall(self,prefix,md,fid) % {{{
    … …  
    5249                        WriteData(fid,prefix,'object',self,'fieldname','spcwatercolumn','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    5350                        WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Double');
    54       outputs = self.requested_outputs;
    55       pos  = find(ismember(outputs,'default'));
    56       if ~isempty(pos),
    57         outputs(pos) = [];  %remove 'default' from outputs
    58         outputs      = [outputs defaultoutputs(self,md)]; %add defaults
    59       end
    60       WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
     51                        outputs = self.requested_outputs;
     52                        pos = find(ismember(outputs,'default'));
     53                        if ~isempty(pos),
     54                                outputs(pos) = []; %remove 'default' from outputs
     55                                outputs      = [outputs defaultoutputs(self,md)]; %add defaults
     56                        end
     57                        WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
     58                end % }}}
     59                function self = extrude(self,md) % {{{
    6160                end % }}}
    6261                function savemodeljs(self,fid,modelname) % {{{

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

    r26059 r26301  
     1import numpy as np
     2
    13from fielddisplay import fielddisplay
    24from checkfield import checkfield
    … …  
    1113    """
    1214
    13     def __init__(self):  # {{{
    14         self.spcwatercolumn = float('NaN')
     15    def __init__(self, *args): #{{{
     16        self.spcwatercolumn = np.nan
    1517        self.stabilization = 0
    1618        self.requested_outputs = []
    1719
    18         self.setdefaultparameters()
     20        nargs = len(args)
     21        if nargs == 0:
     22            self.setdefaultparameters()
     23        elif nargs == 1:
     24            self.setdefaultparameters(args)
     25        else:
     26            raise RuntimeError('constructor not supported')
    1927    #}}}
    2028
    21     def __repr__(self):  # {{{
    22         # TODO:
    23         # - Convert all formatting to calls to <string>.format (see any
    24         #   already converted <class>.__repr__ method for examples)
    25         #
    26         string = '   hydrologyshreve solution parameters:'
    27         string = "%s\n%s" % (string, fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]'))
    28         string = "%s\n%s" % (string, fielddisplay(self, 'stabilization', 'artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.'))
    29         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    30         return string
     29    def __repr__(self): #{{{
     30        s = '   hydrologyshreve solution parameters:\n'
     31        s += '{}\n'.format(fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]'))
     32        s += '{}\n'.format(fielddisplay(self, 'stabilization', 'artificial diffusivity (default: 1). can be more than 1 to increase diffusivity.'))
     33        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
     34        return s
    3135    #}}}
    3236
    33     def extrude(self, md):  # {{{
    34         return self
     37    def setdefaultparameters(self): #{{{
     38        # Type of stabilization to use 0:nothing 1:artificial_diffusivity
     39        self.stabilization = 1
     40        self.requested_outputs = ['default']
    3541    #}}}
    3642
    37     def setdefaultparameters(self):  # {{{
    38         #Type of stabilization to use 0:nothing 1:artificial_diffusivity
    39         self.stabilization = 1
    40         self.requested_outputs = ['default']
    41         return self
    42     #}}}
    43 
    44     def defaultoutputs(self, md):  # {{{
    45         list = ['Watercolumn', 'HydrologyWaterVx', 'HydrologyWaterVy']
    46         return list
    47     #}}}
    48 
    49     def checkconsistency(self, md, solution, analyses):  # {{{
     43    def checkconsistency(self, md, solution, analyses): #{{{
    5044        #Early return
    5145        if 'HydrologyShreveAnalysis' not in analyses or (solution == 'TransientSolution' and not md.transient.ishydrology):
    … …  
    5448        md = checkfield(md, 'fieldname', 'hydrology.spcwatercolumn', 'Inf', 1, 'timeseries', 1)
    5549        md = checkfield(md, 'fieldname', 'hydrology.stabilization', '>=', 0)
    56         md = checkfield(md, 'fieldname', 'hydrology.requested_outputs', 'stringrow', 1)
    5750        return md
    5851    # }}}
    5952
    60     def marshall(self, prefix, md, fid):  # {{{
     53    def defaultoutputs(self, md): #{{{
     54        return ['Watercolumn', 'HydrologyWaterVx', 'HydrologyWaterVy']
     55    #}}}
     56
     57    def marshall(self, prefix, md, fid): #{{{
    6158        WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 2, 'format', 'Integer')
    6259        WriteData(fid, prefix, 'object', self, 'fieldname', 'spcwatercolumn', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    6360        WriteData(fid, prefix, 'object', self, 'fieldname', 'stabilization', 'format', 'Double')
    64         #process requested outputs
    6561        outputs = self.requested_outputs
    6662        indices = [i for i, x in enumerate(outputs) if x == 'default']
    67         if len(indices) > 0:
     63        if len(indices):
    6864            outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
    6965            outputs = outputscopy
    7066        WriteData(fid, prefix, 'data', outputs, 'name', 'md.hydrology.requested_outputs', 'format', 'StringArray')
     67    # }}}
    7168
    72     # }}}
     69    def extrude(self, md): #{{{
     70        return self
     71    #}}}

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

    r26241 r26301  
    2727        end
    2828        methods
    29                 function self = extrude(self,md) % {{{
    30                         self.vx=project3d(md,'vector',self.vx,'type','node');
    31                         self.vy=project3d(md,'vector',self.vy,'type','node');
    32                         self.vz=project3d(md,'vector',self.vz,'type','node');
    33                         self.vel=project3d(md,'vector',self.vel,'type','node');
    34                         self.temperature=project3d(md,'vector',self.temperature,'type','node');
    35                         self.enthalpy=project3d(md,'vector',self.enthalpy,'type','node');
    36                         self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node');
    37                         self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node','layer',1);
    38                         self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1);
    39                         self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1);
    40                         self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1);
    41                         self.sealevel=project3d(md,'vector',self.sealevel,'type','node','layer',1);
    42                         self.bottompressure=project3d(md,'vector',self.bottompressure,'type','node','layer',1);
    43                         self.dsl=project3d(md,'vector',self.dsl,'type','node','layer',1);
    44                         self.str=project3d(md,'vector',self.str,'type','node','layer',1);
    45 
    46                         %Lithostatic pressure by default
    47                         self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z);
    48                 end % }}}
    4929                function self = initialization(varargin) % {{{
    5030                        switch nargin
    … …  
    201181                        end
    202182                end % }}}
     183                function self = extrude(self,md) % {{{
     184                        self.vx=project3d(md,'vector',self.vx,'type','node');
     185                        self.vy=project3d(md,'vector',self.vy,'type','node');
     186                        self.vz=project3d(md,'vector',self.vz,'type','node');
     187                        self.vel=project3d(md,'vector',self.vel,'type','node');
     188                        self.temperature=project3d(md,'vector',self.temperature,'type','node');
     189                        self.enthalpy=project3d(md,'vector',self.enthalpy,'type','node');
     190                        self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node');
     191                        self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node','layer',1);
     192                        self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1);
     193                        self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1);
     194                        self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1);
     195                        self.sealevel=project3d(md,'vector',self.sealevel,'type','node','layer',1);
     196                        self.bottompressure=project3d(md,'vector',self.bottompressure,'type','node','layer',1);
     197                        self.dsl=project3d(md,'vector',self.dsl,'type','node','layer',1);
     198                        self.str=project3d(md,'vector',self.str,'type','node','layer',1);
     199
     200                        %Lithostatic pressure by default
     201                        self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z);
     202                end % }}}
    203203                function savemodeljs(self,fid,modelname) % {{{
    204204

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

    r26241 r26301  
    1111
    1212    Usage:
    13     initialization = initialization()
     13        initialization = initialization()
    1414    """
    1515
    16     def __init__(self):  # {{{
     16    def __init__(self): #{{{
    1717        self.vx = np.nan
    1818        self.vy = np.nan
    1919        self.vz = np.nan
    2020        self.vel = np.nan
    21         self.enthalpy = np.nan
    2221        self.pressure = np.nan
    2322        self.temperature = np.nan
     23        self.enthalpy = np.nan
    2424        self.waterfraction = np.nan
    25         self.watercolumn = np.nan
    2625        self.sediment_head = np.nan
    2726        self.epl_head = np.nan
    2827        self.epl_thickness = np.nan
     28        self.watercolumn = np.nan
    2929        self.hydraulic_potential = np.nan
    3030        self.channelarea = np.nan
    3131        self.sealevel = np.nan
    3232        self.bottompressure = np.nan
     33        self.dsl = np.nan
     34        self.str = np.nan
    3335        self.sample = np.nan
    3436
    35         #set defaults
    3637        self.setdefaultparameters()
    3738    #}}}
    38     def __repr__(self):  # {{{
     39    def __repr__(self): #{{{
    3940        s = '   initial field values:\n'
    4041        s += '{}\n'.format(fielddisplay(self, 'vx', 'x component of velocity [m / yr]'))
    … …  
    5253        s += '{}\n'.format(fielddisplay(self, 'hydraulic_potential', 'Hydraulic potential (for GlaDS) [Pa]'))
    5354        s += '{}\n'.format(fielddisplay(self, 'channelarea', 'subglaciale water channel area (for GlaDS) [m2]'))
    54         s += '{}\n'.format(fielddisplay(self,'sample','Realization of a Gaussian random field'))
     55        s += '{}\n'.format(fielddisplay(self,'sample', 'Realization of a Gaussian random field'))
    5556        return s
    5657    #}}}
    57     def extrude(self, md):  # {{{
    58         self.vx = project3d(md, 'vector', self.vx, 'type', 'node')
    59         self.vy = project3d(md, 'vector', self.vy, 'type', 'node')
    60         self.vz = project3d(md, 'vector', self.vz, 'type', 'node')
    61         self.vel = project3d(md, 'vector', self.vel, 'type', 'node')
    62         self.temperature = project3d(md, 'vector', self.temperature, 'type', 'node')
    63         self.enthalpy = project3d(md, 'vector', self.enthalpy, 'type', 'node')
    64         self.waterfraction = project3d(md, 'vector', self.waterfraction, 'type', 'node')
    65         self.watercolumn = project3d(md, 'vector', self.watercolumn, 'type', 'node')
    66         self.sediment_head = project3d(md, 'vector', self.sediment_head, 'type', 'node', 'layer', 1)
    67         self.epl_head = project3d(md, 'vector', self.epl_head, 'type', 'node', 'layer', 1)
    68         self.epl_thickness = project3d(md, 'vector', self.epl_thickness, 'type', 'node', 'layer', 1)
    69         self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node', 'layer', 1)
    70         self.bottompressure = project3d(md, 'vector', self.bottompressure, 'type', 'node', 'layer', 1)
    71 
    72         #Lithostatic pressure by default
    73         if np.ndim(md.geometry.surface) == 2:
    74             print('Reshaping md.geometry.surface for your convenience but you should fix it in your model set up')
    75             self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface.reshape(-1, ) - md.mesh.z)
    76         else:
    77             self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface - md.mesh.z)
    78 
    79         return self
     58    def setdefaultparameters(self): #{{{
     59        return
    8060    #}}}
    81     def setdefaultparameters(self):  # {{{
    82         return self
    83     #}}}
    84     def checkconsistency(self, md, solution, analyses):  # {{{
     61    def checkconsistency(self, md, solution, analyses): #{{{
    8562        if 'StressbalanceAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.isstressbalance:
    8663            if not np.any(np.logical_or(np.isnan(md.initialization.vx), np.isnan(md.initialization.vy))):
    … …  
    130107        if 'SamplingAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.issampling:
    131108            if np.any(np.isnan(md.initialization.sample)):
    132                 md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
     109                md = checkfield(md, 'fieldname', 'initialization.sample', 'NaN', 1,'Inf', 1, 'size', [md.mesh.numberofvertices])
    133110        return md
    134     # }}}
    135     def marshall(self, prefix, md, fid):  # {{{
     111    #}}}
     112    def marshall(self, prefix, md, fid): #{{{
    136113        yts = md.constants.yts
    137         WriteData(fid, prefix, 'object', self, 'fieldname', 'vx', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
    138         WriteData(fid, prefix, 'object', self, 'fieldname', 'vy', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
    139         WriteData(fid, prefix, 'object', self, 'fieldname', 'vz', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
     114
     115        WriteData(fid, prefix, 'object', self, 'fieldname', 'vx', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1 / yts)
     116        WriteData(fid, prefix, 'object', self, 'fieldname', 'vy', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1 / yts)
     117        WriteData(fid, prefix, 'object', self, 'fieldname', 'vz', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1 / yts)
    140118        WriteData(fid, prefix, 'object', self, 'fieldname', 'pressure', 'format', 'DoubleMat', 'mattype', 1)
    141         WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
     119        WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
    142120        WriteData(fid, prefix, 'object', self, 'fieldname', 'bottompressure', 'format', 'DoubleMat', 'mattype', 1)
     121        WriteData(fid, prefix, 'object', self, 'fieldname', 'str', 'format', 'DoubleMat', 'mattype', 1)
     122        WriteData(fid, prefix, 'object', self, 'fieldname', 'dsl', 'format', 'DoubleMat', 'mattype', 1)
    143123        WriteData(fid, prefix, 'object', self, 'fieldname', 'temperature', 'format', 'DoubleMat', 'mattype', 1)
    144124        WriteData(fid, prefix, 'object', self, 'fieldname', 'waterfraction', 'format', 'DoubleMat', 'mattype', 1)
    … …  
    149129        WriteData(fid, prefix, 'object', self, 'fieldname', 'channelarea', 'format', 'DoubleMat', 'mattype', 1)
    150130        WriteData(fid, prefix, 'object', self, 'fieldname', 'hydraulic_potential', 'format', 'DoubleMat', 'mattype', 1)
    151         WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1)
     131        WriteData(fid, prefix, 'object', self, 'fieldname', 'sample', 'format', 'DoubleMat', 'mattype', 1)
     132
    152133        if md.thermal.isenthalpy:
    153134            if (np.size(self.enthalpy) <= 1):
     135                # Reconstruct enthalpy
    154136                tpmp = md.materials.meltingpoint - md.materials.beta * md.initialization.pressure
    155                 pos = np.nonzero(md.initialization.waterfraction > 0.)[0]
     137                pos = np.where(md.initialization.waterfraction > 0)[0]
    156138                self.enthalpy = md.materials.heatcapacity * (md.initialization.temperature - md.constants.referencetemperature)
    157                 self.enthalpy[pos] = md.materials.heatcapacity * (tpmp[pos].reshape(-1, ) - md.constants.referencetemperature) + md.materials.latentheat * md.initialization.waterfraction[pos].reshape(-1, )
     139                self.enthalpy[pos] = md.materials.heatcapacity * (tpmp[pos].reshape(-1, 1) - md.constants.referencetemperature) + md.materials.latentheat * md.initialization.waterfraction[pos].reshape(-1, 1)
    158140
    159141            WriteData(fid, prefix, 'data', self.enthalpy, 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.initialization.enthalpy')
    160     # }}}
     142    #}}}
     143    def extrude(self, md): #{{{
     144        self.vx = project3d(md, 'vector', self.vx, 'type', 'node')
     145        self.vy = project3d(md, 'vector', self.vy, 'type', 'node')
     146        self.vz = project3d(md, 'vector', self.vz, 'type', 'node')
     147        self.vel = project3d(md, 'vector', self.vel, 'type', 'node')
     148        self.temperature = project3d(md, 'vector', self.temperature, 'type', 'node')
     149        self.enthalpy = project3d(md, 'vector', self.enthalpy, 'type', 'node')
     150        self.waterfraction = project3d(md, 'vector', self.waterfraction, 'type', 'node')
     151        self.watercolumn = project3d(md, 'vector', self.watercolumn, 'type', 'node')
     152        self.sediment_head = project3d(md, 'vector', self.sediment_head, 'type', 'node', 'layer', 1)
     153        self.epl_head = project3d(md, 'vector', self.epl_head, 'type', 'node', 'layer', 1)
     154        self.epl_thickness = project3d(md, 'vector', self.epl_thickness, 'type', 'node', 'layer', 1)
     155        self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node', 'layer', 1)
     156        self.bottompressure = project3d(md, 'vector', self.bottompressure, 'type', 'node', 'layer', 1)
     157
     158        # Lithostatic pressure by default
     159        if np.ndim(md.geometry.surface) == 2:
     160            print('Reshaping md.geometry.surface for your convenience but you should fix it in your model set up')
     161            self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface.reshape(-1, 1) - md.mesh.z)
     162        else:
     163            self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface - md.mesh.z)
     164
     165        return self
     166    #}}}

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

    r26272 r26301  
    1616               
    1717                %tidal love numbers for computing rotational feedback:
    18                 th          = []; 
    19                 tk          = []; 
    20                 tl          = []; 
    21                 tk2secular  = 0;  %deg 2 secular number.
     18                th          = [];
     19                tk          = [];
     20                tl          = [];
     21                tk2secular  = 0; %deg 2 secular number.
    2222
    2323                %time/frequency for visco-elastic love numbers
    2424                timefreq    = [];
    25                 istime      = 1; 
     25                istime      = 1;
    2626
    2727        end
    … …  
    3333                        self=setdefaultparameters(self,maxdeg,referenceframe);
    3434                end % }}}
     35                function disp(self) % {{{
     36                        disp(sprintf('   lovenumbers parameters:'));
     37
     38                        fielddisplay(self,'h','load Love number for radial displacement');
     39                        fielddisplay(self,'k','load Love number for gravitational potential perturbation');
     40                        fielddisplay(self,'l','load Love number for horizontal displacements');
     41
     42                        fielddisplay(self,'th','tidal load Love number (deg 2)');
     43                        fielddisplay(self,'tk','tidal load Love number (deg 2)');
     44                        fielddisplay(self,'tl','tidal load Love number (deg 2)');
     45                        fielddisplay(self,'tk2secular','secular fluid Love number');
     46
     47                        fielddisplay(self,'istime','time (default: 1) or frequency love numbers (0)');
     48                        fielddisplay(self,'timefreq','time/frequency vector (yr or 1/yr)');
     49
     50                end % }}}
    3551                function self = setdefaultparameters(self,maxdeg,referenceframe) % {{{
    36                
     52
    3753                        %initialize love numbers:
    3854                        self.h=getlovenumbers('type','loadingverticaldisplacement','referenceframe',referenceframe,'maxdeg',maxdeg);
    … …  
    4561                        %secular fluid love number:
    4662                        self.tk2secular=0.942;
    47                        
     63
    4864                        %time:
    4965                        self.istime=1; %temporal love numbers by default
    … …  
    7389                                error('lovenumbers error message: love numbers should be provided at the same level of accuracy');
    7490                        end
    75                        
     91
    7692                        ntf=length(self.timefreq);
    7793                        if( size(self.h,2) ~= ntf | size(self.k,2) ~= ntf | size(self.l,2) ~= ntf | size(self.th,2) ~= ntf | size(self.tk,2) ~= ntf | size(self.tl,2) ~= ntf ),
    … …  
    8298                function list=defaultoutputs(self,md) % {{{
    8399                        list = {};
    84                 end % }}}
    85                 function disp(self) % {{{
    86                         disp(sprintf('   lovenumbers parameters:'));
    87 
    88                         fielddisplay(self,'h','load Love number for radial displacement');
    89                         fielddisplay(self,'k','load Love number for gravitational potential perturbation');
    90                         fielddisplay(self,'l','load Love number for horizontal displacements');
    91 
    92                         fielddisplay(self,'th','tidal load Love number (deg 2)');
    93                         fielddisplay(self,'tk','tidal load Love number (deg 2)');
    94                         fielddisplay(self,'tl','tidal load Love number (deg 2)');
    95                         fielddisplay(self,'tk2secular','secular fluid Love number');
    96                        
    97                         fielddisplay(self,'istime','time (default=1) or frequency love numbers (0)');
    98                         fielddisplay(self,'timefreq','time/frequency vector (yr or 1/yr)');
    99 
    100100                end % }}}
    101101                function marshall(self,prefix,md,fid) % {{{

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

    r26178 r26301  
    1010
    1111    Usage:
    12         lovenumbers = lovenumbers()  #will setup love numbers deg 1001 by default
    13         lovenumbers = lovenumbers('maxdeg', 10001);   #supply numbers of degrees required (here, 10001)
     12        lovenumbers = lovenumbers()
     13        lovenumbers = lovenumbers('maxdeg', 10000, 'referenceframe', 'CF');
     14
     15    Choose numbers of degrees required (1000 by default) and reference frame
     16    (between CF and CM; CM by default)
    1417    """
    1518
    … …  
    2629        self.tk2secular = 0  # deg 2 secular number
    2730
     31        # Time/frequency for visco-elastic love numbers
     32        self.timefreq = []
     33        self.istime = 1
     34
    2835        options = pairoptions(*args)
    2936        maxdeg = options.getfieldvalue('maxdeg', 1000)
    3037        referenceframe = options.getfieldvalue('referenceframe', 'CM')
    3138        self.setdefaultparameters(maxdeg, referenceframe)
     39    #}}}
     40
     41    def __repr__(self):  #{{{
     42        s = '   lovenumbers parameters:\n'
     43        s += '{}\n'.format(fielddisplay(self, 'h', 'load Love number for radial displacement'))
     44        s += '{}\n'.format(fielddisplay(self, 'k', 'load Love number for gravitational potential perturbation'))
     45        s += '{}\n'.format(fielddisplay(self, 'l', 'load Love number for horizontal displacements'))
     46        s += '{}\n'.format(fielddisplay(self, 'th', 'tidal load Love number (deg 2)'))
     47        s += '{}\n'.format(fielddisplay(self, 'tk', 'tidal load Love number (deg 2)'))
     48        s += '{}\n'.format(fielddisplay(self, 'tl', 'tidal load Love number (deg 2)'))
     49        s += '{}\n'.format(fielddisplay(self, 'tk2secular', 'secular fluid Love number'))
     50        s += '{}\n'.format(fielddisplay(self, 'istime', 'time (default: 1) or frequency love numbers (0)'))
     51        s += '{}\n'.format(fielddisplay(self, 'timefreq', 'time/frequency vector (yr or 1/yr)'))
     52        return s
    3253    #}}}
    3354
    … …  
    4364        # Secular fluid love number
    4465        self.tk2secular = 0.942
     66
     67        # Time
     68        self.istime = 1 # Temporal love numbers by default
     69        self.timefreq = 0 # Elastic case by default
    4570        return self
    4671    #}}}
    4772
    4873    def checkconsistency(self, md, solution, analyses):  #{{{
    49         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
     74        if ('SealevelchangeAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
    5075            return
    5176        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.h', 'NaN', 1, 'Inf', 1)
    … …  
    5681        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.tk', 'NaN', 1, 'Inf', 1)
    5782        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.tk2secular', 'NaN', 1, 'Inf', 1)
     83        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.timefreq', 'NaN', 1, 'Inf', 1)
     84        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.istime', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
     85
    5886        # Check that love numbers are provided at the same level of accuracy
    5987        if (self.h.shape[0] != self.k.shape[0]) or (self.h.shape[0] != self.l.shape[0]):
    6088            raise ValueError('lovenumbers error message: love numbers should be provided at the same level of accuracy')
     89
     90        ntf = len(self.timefreq)
     91        if (self.h.shape[1] != ntf or self.k.shape[1] != ntf or self.l.shape[1] != ntf or self.th.shape[1] != ntf or self.tk.shape[1] != ntf or self.tl.shape[1] != ntf):
     92            raise ValueError('lovenumbers error message: love numbers should have as many time/frequency steps as the time/frequency vector')
     93
    6194        return md
    6295    #}}}
    … …  
    6497    def defaultoutputs(self, md):  #{{{
    6598        return[]
    66     #}}}
    67 
    68     def __repr__(self):  #{{{
    69         s = '   lovenumbers parameters:\n'
    70         s += '{}\n'.format(fielddisplay(self, 'h', 'load Love number for radial displacement'))
    71         s += '{}\n'.format(fielddisplay(self, 'k', 'load Love number for gravitational potential perturbation'))
    72         s += '{}\n'.format(fielddisplay(self, 'l', 'load Love number for horizontal displacements'))
    73         s += '{}\n'.format(fielddisplay(self, 'th', 'tidal load Love number (deg 2)'))
    74         s += '{}\n'.format(fielddisplay(self, 'tk', 'tidal load Love number (deg 2)'))
    75         s += '{}\n'.format(fielddisplay(self, 'tk2secular', 'secular fluid Love number'))
    76         return s
    7799    #}}}
    78100
    … …  
    86108        WriteData(fid, prefix, 'object', self, 'fieldname', 'tl', 'name', 'md.solidearth.lovenumbers.tl', 'format', 'DoubleMat', 'mattype', 1)
    87109        WriteData(fid, prefix, 'object', self, 'data', self.tk2secular, 'fieldname', 'lovenumbers.tk2secular', 'format', 'Double')
     110
     111        if (self.istime):
     112            scale = md.constants.yts
     113        else:
     114            scale = 1.0 / md.constants.yts
     115        WriteData(fid, prefix, 'object', self, 'fieldname', 'istime', 'name', 'md.solidearth.lovenumbers.istime', 'format', 'Boolean')
     116        WriteData(fid, prefix, 'object', self, 'fieldname', 'timefreq', 'name', 'md.solidearth.lovenumbers.timefreq', 'format', 'DoubleMat', 'mattype', 1, 'scale', scale);
    88117    #}}}
    89118

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

    r26299 r26301  
    204204                                        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)');
    205205                                        fielddisplay(self,'ebm_taul','array describing each layer''s starting period for transient relaxation, only for EBM rheology  (numlayers) [s]');
    206                                         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]');
     206                                        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]');
    207207
    208208
    … …  
    256256                                                        error('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with rheologymodel choice');
    257257                                                end
    258                                                 if md.materials.rheologymodel(i)==2 & (isnan(md.materials.ebm_alpha(i)) | isnan(md.materials.ebm_delta(i)) | isnan(md.materials.ebm_taul(i)) | isnan(md.materials.ebm_tauh(i)) ),
     258                                                if md.materials.rheologymodel(i)==2 & (isnan(md.materials.ebm_alpha(i)) | isnan(md.materials.ebm_delta(i)) | isnan(md.materials.ebm_taul(i)) | isnan(md.materials.ebm_tauh(i))),
    259259                                                        error('materials checkconsistency error message: Litho ebm_alpha, ebm_delta, ebm_taul or ebm_tauh has NaN values, inconsistent with rheologymodel choice');
    260260                                                end
    … …  
    284284                        %1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum
    285285                        WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(self.nature),'format','IntMat','mattype',3);
    286                         WriteData(fid,prefix,'name','md.materials.type','data',5,'format','Integer'); %DANGER, this can evolve if you have classes.
     286                        WriteData(fid,prefix,'name','md.materials.type','data',5,'format','Integer'); %DANGER: this can evolve if you have classes.
    287287                        for i=1:length(self.nature),
    288288                                nat=self.nature{i};
    … …  
    326326                                        end
    327327                                        earth_density=earth_density/self.radius(self.numlayers+1)^3;
     328                                        disp(earth_density)
    328329                                        self.earth_density=earth_density;
    329330                                case 'hydro'
    … …  
    541542                                vs =  vs + t1*(r2^3) + t2*(r2^4) + t3*(r2^5) + t4*(r2^6) - ( t1*(r1^3) + t2*(r1^4) + t3*(r1^5) + t4*(r1^6) );
    542543
    543                         end 
     544                        end
    544545                        ro = ro*3 / (rad(j+1)^3-rad(j)^3);
    545546                        vp = vp*3 /(rad(j+1)^3-rad(j)^3);
    … …  
    557558
    558559                end % }}}
    559    
    560560        end
    561561end
    … …  
    584584        end
    585585end % }}}
    586 
    587    

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

    r26299 r26301  
    2525                raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
    2626
     27        # Start filling in the dynamic fields (not truly dynamic under Python)
    2728        for i in range(len(self.nature)):
    2829            nat = self.nature[i]
    … …  
    6768        self.earth_density = 0
    6869
    69         # Set default parameters
    7070        self.setdefaultparameters()
    7171    #}}}
    … …  
    105105                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)'))
    106106                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]'))
     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]'))
    108108
    109109                s += '{}\n'.format(fielddisplay(self, 'rheologymodel', 'array describing whether we adopt a Maxwell (0), Burgers (1) or EBM (2) rheology (default: 0)'))
    110110                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)'))
     111                s += '{}\n'.format(fielddisplay(self, 'issolid', 'array describing whether the layer is solid or liquid (default: 1) (numlayers)'))
    112112            elif nat == 'hydro':
    113113                s += 'Hydro:\n'
    … …  
    127127            if nat == 'ice':
    128128                # Ice density (kg/m^3)
    129                 self.rho_ice = 917.
     129                self.rho_ice = 917
    130130
    131131                # Ocean water density (kg/m^3)
    132                 self.rho_water = 1023.
     132                self.rho_water = 1023
    133133
    134134                # Fresh water density (kg/m^3)
    135                 self.rho_freshwater = 1000.
     135                self.rho_freshwater = 1000
    136136
    137137                # Water viscosity (N.s/m^2)
    … …  
    139139
    140140                # Ice heat capacity cp (J/kg/K)
    141                 self.heatcapacity = 2093.
     141                self.heatcapacity = 2093
    142142
    143143                # Ice latent heat of fusion L (J/kg)
    … …  
    160160
    161161                # Mixed layer (ice-water interface) heat capacity (J/kg/K)
    162                 self.mixed_layer_capacity = 3974.
     162                self.mixed_layer_capacity = 3974
    163163
    164164                # Thermal exchange velocity (ice-water interface) (m/s)
    … …  
    193193                self.ebm_tauh = [np.nan, np.nan]
    194194                self.rheologymodel = [0, 0]
    195                 self.density = [5.51e3, 5.50e3]  # (Pa) # Mantle and lithosphere density [kg/m^3]
    196                 self.issolid = [1, 1]  # Is layer solid or liquid?
     195                self.density = [5.51e3, 5.50e3] # (Pa) # Mantle and lithosphere density [kg/m^3]
     196                self.issolid = [1, 1] # Is layer solid or liquid?
    197197            elif nat == 'hydro':
    198198                # Ice density (kg/m^3)
    199                 self.rho_ice = 917.
     199                self.rho_ice = 917
    200200
    201201                # Ocean water density (kg/m^3)
    202                 self.rho_water = 1023.
     202                self.rho_water = 1023
    203203
    204204                # Fresh water density (kg/m^3)
    205                 self.rho_freshwater = 1000.
     205                self.rho_freshwater = 1000
    206206            else:
    207207                raise RuntimeError("materials setdefaultparameters error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
    … …  
    251251                ind = np.where(md.materials.issolid == 0)[0]
    252252                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 #{} but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.'.format(i))
     253                    raise RuntimeError('Fluid layers detected at layers #' + indices + ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.')
    254254
    255255            elif nat == 'hydro':
    … …  
    267267        #1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum
    268268        WriteData(fid, prefix, 'name', 'md.materials.nature', 'data', naturetointeger(self.nature), 'format', 'IntMat', 'mattype', 3)
    269         WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER, this can evolve if you have classes
     269        WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER: this can evolve if you have classes
    270270        for i in range(len(self.nature)):
    271271            nat = self.nature[i]
    … …  
    307307                    earth_density = earth_density + (pow(self.radius[i + 1], 3) - pow(self.radius[i], 3)) * self.density[i]
    308308                earth_density = earth_density / pow(self.radius[self.numlayers], 3)
     309                print(earth_density)
    309310                self.earth_density = earth_density
    310311            elif nat == 'hydro':

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

    r26213 r26301  
    3939                steadystate      = 0;
    4040                transient        = 0;
    41                 levelset             = 0;
     41                levelset         = 0;
    4242                calving          = 0;
    4343                frontalforcings  = 0;
    44                 love                    = 0;
     44                love            = 0;
    4545                esa              = 0;
    4646                sampling         = 0;
    … …  
    4949                inversion        = 0;
    5050                qmu              = 0;
    51                 amr                              = 0;
     51                amr              = 0;
    5252                results          = 0;
    5353                outputdefinition = 0;
    … …  
    207207                end
    208208                %}}}
     209                function disp(self) % {{{
     210                        disp(sprintf('%19s: %-22s -- %s','mesh'            ,['[1x1 ' class(self.mesh) ']'],'mesh properties'));
     211                        disp(sprintf('%19s: %-22s -- %s','mask'            ,['[1x1 ' class(self.mask) ']'],'defines grounded and floating elements'));
     212                        disp(sprintf('%19s: %-22s -- %s','geometry'        ,['[1x1 ' class(self.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...'));
     213                        disp(sprintf('%19s: %-22s -- %s','constants'       ,['[1x1 ' class(self.constants) ']'],'physical constants'));
     214                        disp(sprintf('%19s: %-22s -- %s','smb'             ,['[1x1 ' class(self.smb) ']'],'surface mass balance'));
     215                        disp(sprintf('%19s: %-22s -- %s','basalforcings'   ,['[1x1 ' class(self.basalforcings) ']'],'bed forcings'));
     216                        disp(sprintf('%19s: %-22s -- %s','materials'       ,['[1x1 ' class(self.materials) ']'],'material properties'));
     217                        disp(sprintf('%19s: %-22s -- %s','damage'          ,['[1x1 ' class(self.damage) ']'],'parameters for damage evolution solution'));
     218                        disp(sprintf('%19s: %-22s -- %s','friction'        ,['[1x1 ' class(self.friction) ']'],'basal friction/drag properties'));
     219                        disp(sprintf('%19s: %-22s -- %s','flowequation'    ,['[1x1 ' class(self.flowequation) ']'],'flow equations'));
     220                        disp(sprintf('%19s: %-22s -- %s','timestepping'    ,['[1x1 ' class(self.timestepping) ']'],'time stepping for transient models'));
     221                        disp(sprintf('%19s: %-22s -- %s','initialization'  ,['[1x1 ' class(self.initialization) ']'],'initial guess/state'));
     222                        disp(sprintf('%19s: %-22s -- %s','rifts'           ,['[1x1 ' class(self.rifts) ']'],'rifts properties'));
     223                        disp(sprintf('%19s: %-22s -- %s','solidearth'      ,['[1x1 ' class(self.solidearth) ']'],'solidearth inputs and settings'));
     224                        disp(sprintf('%19s: %-22s -- %s','dsl'             ,['[1x1 ' class(self.dsl) ']'],'dynamic sea-level '));
     225                        disp(sprintf('%19s: %-22s -- %s','debug'           ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)'));
     226                        disp(sprintf('%19s: %-22s -- %s','verbose'         ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve'));
     227                        disp(sprintf('%19s: %-22s -- %s','settings'        ,['[1x1 ' class(self.settings) ']'],'settings properties'));
     228                        disp(sprintf('%19s: %-22s -- %s','toolkits'        ,['[1x1 ' class(self.toolkits) ']'],'PETSc options for each solution'));
     229                        disp(sprintf('%19s: %-22s -- %s','cluster'         ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of CPUs...)'));
     230                        disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(self.balancethickness) ']'],'parameters for balancethickness solution'));
     231                        disp(sprintf('%19s: %-22s -- %s','stressbalance'   ,['[1x1 ' class(self.stressbalance) ']'],'parameters for stressbalance solution'));
     232                        disp(sprintf('%19s: %-22s -- %s','groundingline'   ,['[1x1 ' class(self.groundingline) ']'],'parameters for groundingline solution'));
     233                        disp(sprintf('%19s: %-22s -- %s','hydrology'       ,['[1x1 ' class(self.hydrology) ']'],'parameters for hydrology solution'));
     234                        disp(sprintf('%19s: %-22s -- %s','masstransport'   ,['[1x1 ' class(self.masstransport) ']'],'parameters for masstransport solution'));
     235                        disp(sprintf('%19s: %-22s -- %s','thermal'         ,['[1x1 ' class(self.thermal) ']'],'parameters for thermal solution'));
     236                        disp(sprintf('%19s: %-22s -- %s','steadystate'     ,['[1x1 ' class(self.steadystate) ']'],'parameters for steadystate solution'));
     237                        disp(sprintf('%19s: %-22s -- %s','transient'       ,['[1x1 ' class(self.transient) ']'],'parHwoameters for transient solution'));
     238                        disp(sprintf('%19s: %-22s -- %s','levelset'        ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)'));
     239                        disp(sprintf('%19s: %-22s -- %s','calving'         ,['[1x1 ' class(self.calving) ']'],'parameters for calving'));
     240                        disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings'));
     241                        disp(sprintf('%19s: %-22s -- %s','esa'             ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution'));
     242                        disp(sprintf('%19s: %-22s -- %s','love'            ,['[1x1 ' class(self.love) ']'],'parameters for love solution'));
     243                        disp(sprintf('%19s: %-22s -- %s','sampling'        ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler'));
     244                        disp(sprintf('%19s: %-22s -- %s','autodiff'        ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters'));
     245                        disp(sprintf('%19s: %-22s -- %s','inversion'       ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods'));
     246                        disp(sprintf('%19s: %-22s -- %s','qmu'             ,['[1x1 ' class(self.qmu) ']'],'Dakota properties'));
     247                        disp(sprintf('%19s: %-22s -- %s','amr'             ,['[1x1 ' class(self.amr) ']'],'adaptive mesh refinement properties'));
     248                        disp(sprintf('%19s: %-22s -- %s','outputdefinition',['[1x1 ' class(self.outputdefinition) ']'],'output definition'));
     249                        disp(sprintf('%19s: %-22s -- %s','results'         ,['[1x1 ' class(self.results) ']'],'model results'));
     250                        disp(sprintf('%19s: %-22s -- %s','radaroverlay'    ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay'));
     251                        disp(sprintf('%19s: %-22s -- %s','miscellaneous'   ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields'));
     252                end % }}}
    209253                function md = setdefaultparameters(md,planet) % {{{
    210254
    … …  
    11331177                                md.mesh.average_vertex_connectivity=100;
    11341178                        end
    1135                         end % }}}
     1179                end % }}}
    11361180                function md = structtomodel(md,structmd) % {{{
    11371181
    … …  
    15841628                        end
    15851629                end % }}}
    1586                 function disp(self) % {{{
    1587                         disp(sprintf('%19s: %-22s -- %s','mesh'            ,['[1x1 ' class(self.mesh) ']'],'mesh properties'));
    1588                         disp(sprintf('%19s: %-22s -- %s','mask'            ,['[1x1 ' class(self.mask) ']'],'defines grounded and floating elements'));
    1589                         disp(sprintf('%19s: %-22s -- %s','geometry'        ,['[1x1 ' class(self.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...'));
    1590                         disp(sprintf('%19s: %-22s -- %s','constants'       ,['[1x1 ' class(self.constants) ']'],'physical constants'));
    1591                         disp(sprintf('%19s: %-22s -- %s','smb'             ,['[1x1 ' class(self.smb) ']'],'surface mass balance'));
    1592                         disp(sprintf('%19s: %-22s -- %s','basalforcings'   ,['[1x1 ' class(self.basalforcings) ']'],'bed forcings'));
    1593                         disp(sprintf('%19s: %-22s -- %s','materials'       ,['[1x1 ' class(self.materials) ']'],'material properties'));
    1594                         disp(sprintf('%19s: %-22s -- %s','damage'          ,['[1x1 ' class(self.damage) ']'],'parameters for damage evolution solution'));
    1595                         disp(sprintf('%19s: %-22s -- %s','friction'        ,['[1x1 ' class(self.friction) ']'],'basal friction/drag properties'));
    1596                         disp(sprintf('%19s: %-22s -- %s','flowequation'    ,['[1x1 ' class(self.flowequation) ']'],'flow equations'));
    1597                         disp(sprintf('%19s: %-22s -- %s','timestepping'    ,['[1x1 ' class(self.timestepping) ']'],'time stepping for transient models'));
    1598                         disp(sprintf('%19s: %-22s -- %s','initialization'  ,['[1x1 ' class(self.initialization) ']'],'initial guess/state'));
    1599                         disp(sprintf('%19s: %-22s -- %s','rifts'           ,['[1x1 ' class(self.rifts) ']'],'rifts properties'));
    1600                         disp(sprintf('%19s: %-22s -- %s','solidearth'      ,['[1x1 ' class(self.solidearth) ']'],'solidearth inputs and settings'));
    1601                         disp(sprintf('%19s: %-22s -- %s','dsl'             ,['[1x1 ' class(self.dsl) ']'],'dynamic sea-level '));
    1602                         disp(sprintf('%19s: %-22s -- %s','debug'           ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)'));
    1603                         disp(sprintf('%19s: %-22s -- %s','verbose'         ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve'));
    1604                         disp(sprintf('%19s: %-22s -- %s','settings'        ,['[1x1 ' class(self.settings) ']'],'settings properties'));
    1605                         disp(sprintf('%19s: %-22s -- %s','toolkits'        ,['[1x1 ' class(self.toolkits) ']'],'PETSc options for each solution'));
    1606                         disp(sprintf('%19s: %-22s -- %s','cluster'         ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of CPUs...)'));
    1607                         disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(self.balancethickness) ']'],'parameters for balancethickness solution'));
    1608                         disp(sprintf('%19s: %-22s -- %s','stressbalance'   ,['[1x1 ' class(self.stressbalance) ']'],'parameters for stressbalance solution'));
    1609                         disp(sprintf('%19s: %-22s -- %s','groundingline'   ,['[1x1 ' class(self.groundingline) ']'],'parameters for groundingline solution'));
    1610                         disp(sprintf('%19s: %-22s -- %s','hydrology'       ,['[1x1 ' class(self.hydrology) ']'],'parameters for hydrology solution'));
    1611                         disp(sprintf('%19s: %-22s -- %s','masstransport'   ,['[1x1 ' class(self.masstransport) ']'],'parameters for masstransport solution'));
    1612                         disp(sprintf('%19s: %-22s -- %s','thermal'         ,['[1x1 ' class(self.thermal) ']'],'parameters for thermal solution'));
    1613                         disp(sprintf('%19s: %-22s -- %s','steadystate'     ,['[1x1 ' class(self.steadystate) ']'],'parameters for steadystate solution'));
    1614                         disp(sprintf('%19s: %-22s -- %s','transient'       ,['[1x1 ' class(self.transient) ']'],'parHwoameters for transient solution'));
    1615                         disp(sprintf('%19s: %-22s -- %s','levelset'        ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)'));
    1616                         disp(sprintf('%19s: %-22s -- %s','calving'         ,['[1x1 ' class(self.calving) ']'],'parameters for calving'));
    1617                         disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings'));
    1618                         disp(sprintf('%19s: %-22s -- %s','esa'             ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution'));
    1619                         disp(sprintf('%19s: %-22s -- %s','love'            ,['[1x1 ' class(self.love) ']'],'parameters for love solution'));
    1620                         disp(sprintf('%19s: %-22s -- %s','sampling'        ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler'));
    1621                         disp(sprintf('%19s: %-22s -- %s','autodiff'        ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters'));
    1622                         disp(sprintf('%19s: %-22s -- %s','inversion'       ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods'));
    1623                         disp(sprintf('%19s: %-22s -- %s','qmu'             ,['[1x1 ' class(self.qmu) ']'],'Dakota properties'));
    1624                         disp(sprintf('%19s: %-22s -- %s','amr'             ,['[1x1 ' class(self.amr) ']'],'adaptive mesh refinement properties'));
    1625                         disp(sprintf('%19s: %-22s -- %s','outputdefinition',['[1x1 ' class(self.outputdefinition) ']'],'output definition'));
    1626                         disp(sprintf('%19s: %-22s -- %s','results'         ,['[1x1 ' class(self.results) ']'],'model results'));
    1627                         disp(sprintf('%19s: %-22s -- %s','radaroverlay'    ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay'));
    1628                         disp(sprintf('%19s: %-22s -- %s','miscellaneous'   ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields'));
    1629                 end % }}}
    16301630                function memory(self) % {{{
    16311631
    1632                 disp(sprintf('\nMemory imprint:\n'));
    1633 
    1634                 fields=properties('model');
    1635                 mem=0;
    1636 
    1637                 for i=1:length(fields),
    1638                         field=self.(fields{i});
    1639                         s=whos('field');
    1640                         mem=mem+s.bytes/1e6;
    1641                         disp(sprintf('%19s: %6.2f Mb',fields{i},s.bytes/1e6));
     1632                        disp(sprintf('\nMemory imprint:\n'));
     1633
     1634                        fields=properties('model');
     1635                        mem=0;
     1636
     1637                        for i=1:length(fields),
     1638                                field=self.(fields{i});
     1639                                s=whos('field');
     1640                                mem=mem+s.bytes/1e6;
     1641                                disp(sprintf('%19s: %6.2f Mb',fields{i},s.bytes/1e6));
     1642                        end
     1643                        disp(sprintf('%19s--%10s','--------------','--------------'));
     1644                        disp(sprintf('%19s: %g Mb','Total',mem));
    16421645                end
    1643                 disp(sprintf('%19s--%10s','--------------','--------------'));
    1644                 disp(sprintf('%19s: %g Mb','Total',mem));
    1645                 end % }}}
     1646                % }}}
    16461647                function netcdf(self,filename) % {{{
    1647                 %NETCDF - save model as netcdf
    1648                 %
    1649                 %   Usage:
    1650                 %      netcdf(md,filename)
    1651                 %
    1652                 %   Example:
    1653                 %      netcdf(md,'model.nc');
    1654 
    1655                 disp('Saving model as NetCDF');
    1656                 %1. Create NetCDF file
    1657                 ncid=netcdf.create(filename,'CLOBBER');
    1658                 netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.4');
    1659                 netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Title',['ISSM model (' self.miscellaneous.name ')']);
    1660                 netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Author',getenv('USER'));
    1661                 netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Date',datestr(now));
    1662 
    1663                 %Preallocate variable id, needed to write variables in netcdf file
    1664                 var_id=zeros(1000,1);%preallocate
    1665 
    1666                 for step=1:2,
    1667                         counter=0;
    1668                         [var_id,counter]=structtonc(ncid,'md',self,0,var_id,counter,step);
    1669                         if step==1, netcdf.endDef(ncid); end
    1670                 end
    1671 
    1672                 if counter>1000,
    1673                         warning(['preallocation of var_id need to be updated from ' num2str(1000) ' to ' num2str(counter)]);
    1674                 end
    1675 
    1676                 netcdf.close(ncid)
     1648                        %NETCDF - save model as netcdf
     1649                        %
     1650                        %   Usage:
     1651                        %      netcdf(md,filename)
     1652                        %
     1653                        %   Example:
     1654                        %      netcdf(md,'model.nc');
     1655
     1656                        disp('Saving model as NetCDF');
     1657                        %1. Create NetCDF file
     1658                        ncid=netcdf.create(filename,'CLOBBER');
     1659                        netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.4');
     1660                        netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Title',['ISSM model (' self.miscellaneous.name ')']);
     1661                        netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Author',getenv('USER'));
     1662                        netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Date',datestr(now));
     1663
     1664                        %Preallocate variable id, needed to write variables in netcdf file
     1665                        var_id=zeros(1000,1);%preallocate
     1666
     1667                        for step=1:2,
     1668                                counter=0;
     1669                                [var_id,counter]=structtonc(ncid,'md',self,0,var_id,counter,step);
     1670                                if step==1, netcdf.endDef(ncid); end
     1671                        end
     1672
     1673                        if counter>1000,
     1674                                warning(['preallocation of var_id need to be updated from ' num2str(1000) ' to ' num2str(counter)]);
     1675                        end
     1676
     1677                        netcdf.close(ncid)
    16771678                end % }}}
    16781679                function xylim(self) % {{{
    … …  
    16821683                end % }}}
    16831684                function md=upload(md) % {{{
    1684                 %the goal of this routine is to upload the model onto a server, and to empty it.
    1685                 %So first, save the model with a unique name and upload the file to the server:
    1686                 random_part=fix(rand(1)*10000);
    1687                 id=[md.miscellaneous.name '-' regexprep(datestr(now),'[^\w'']','') '-' num2str(random_part)  '-' getenv('USER') '-' oshostname() '.upload'];
    1688                 eval(['save ' id ' md']);
    1689 
    1690                 %Now, upload the file:
    1691                 issmscpout(md.settings.upload_server,md.settings.upload_path,md.settings.upload_login,md.settings.upload_port,{id},1);
    1692 
    1693                 %Now, empty this model of everything except settings, and record name of file we just uploaded!
    1694                 settings_back=md.settings;
    1695                 md=model();
    1696                 md.settings=settings_back;
    1697                 md.settings.upload_filename=id;
    1698 
    1699                 %get locally rid of file that was uploaded
    1700                 eval(['delete ' id]);
     1685                        %the goal of this routine is to upload the model onto a server, and to empty it.
     1686                        %So first, save the model with a unique name and upload the file to the server:
     1687                        random_part=fix(rand(1)*10000);
     1688                        id=[md.miscellaneous.name '-' regexprep(datestr(now),'[^\w'']','') '-' num2str(random_part)  '-' getenv('USER') '-' oshostname() '.upload'];
     1689                        eval(['save ' id ' md']);
     1690
     1691                        %Now, upload the file:
     1692                        issmscpout(md.settings.upload_server,md.settings.upload_path,md.settings.upload_login,md.settings.upload_port,{id},1);
     1693
     1694                        %Now, empty this model of everything except settings, and record name of file we just uploaded!
     1695                        settings_back=md.settings;
     1696                        md=model();
     1697                        md.settings=settings_back;
     1698                        md.settings.upload_filename=id;
     1699
     1700                        %get locally rid of file that was uploaded
     1701                        eval(['delete ' id]);
    17011702
    17021703                end % }}}
    17031704                function md=download(md) % {{{
    17041705
    1705                 %the goal of this routine is to download the internals of the current model from a server, because
    1706                 %this model is empty, except for the settings which tell us where to go and find this model!
    1707 
    1708                 %Download the file:
    1709                 issmscpin(md.settings.upload_server, md.settings.upload_login, md.settings.upload_port, md.settings.upload_path, {md.settings.upload_filename});
    1710 
    1711                 name=md.settings.upload_filename;
    1712 
    1713                 %Now, load this model:
    1714                 md=loadmodel(md.settings.upload_filename);
    1715 
    1716                 %get locally rid of file that was downloaded
    1717                 eval(['delete ' name]);
     1706                        %the goal of this routine is to download the internals of the current model from a server, because
     1707                        %this model is empty, except for the settings which tell us where to go and find this model!
     1708
     1709                        %Download the file:
     1710                        issmscpin(md.settings.upload_server, md.settings.upload_login, md.settings.upload_port, md.settings.upload_path, {md.settings.upload_filename});
     1711
     1712                        name=md.settings.upload_filename;
     1713
     1714                        %Now, load this model:
     1715                        md=loadmodel(md.settings.upload_filename);
     1716
     1717                        %get locally rid of file that was downloaded
     1718                        eval(['delete ' name]);
    17181719
    17191720                end % }}}

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

    r26185 r26301  
    7676
    7777class model(object):
    78     #properties
     78    """MODEL - class definition
     79
     80    Usage:
     81        md = model()
     82    """
     83
    7984    def __init__(self, *args): #{{{
    8085        self.mesh = None
    8186        self.mask = None
     87
     88        self.geometry = None
    8289        self.constants = None
    83         self.geometry = None
    84         self.initialization = None
    8590        self.smb = None
    8691        self.basalforcings = None
    87         self.friction = None
    88         self.rifts = None
    89         self.solidearth = None
    90         self.dsl = None
    91         self.timestepping = None
    92         self.groundingline = None
    9392        self.materials = None
    9493        self.damage = None
     94        self.friction = None
    9595        self.flowequation = None
     96        self.timestepping = None
     97        self.initialization = None
     98        self.rifts = None
     99        self.dsl = None
     100        self.solidearth = None
     101
    96102        self.debug = None
    97103        self.verbose = None
    … …  
    99105        self.toolkits = None
    100106        self.cluster = None
     107
    101108        self.balancethickness = None
    102109        self.stressbalance = None
     110        self.groundingline = None
    103111        self.hydrology = None
    104112        self.masstransport = None
    … …  
    112120        self.esa = None
    113121        self.sampling = None
     122
    114123        self.autodiff = None
    115124        self.inversion = None
    116125        self.qmu = None
    117126        self.amr = None
    118         self.radaroverlay = None
    119127        self.results = None
    120128        self.outputdefinition = None
     129        self.radaroverlay = None
    121130        self.miscellaneous = None
    122131        self.private = None
    123132
    124         nargs = len(args)
    125 
    126         if nargs == 0:
     133        if len(args) == 0:
    127134            self.setdefaultparameters('earth')
    128135        else:
    129             self.setdefaultparameters(args[0])
    130136            options = pairoptions(*args)
    131137            planet = options.getfieldvalue('planet', 'earth')
    132138            self.setdefaultparameters(planet)
    133 #}}}
    134 
    135     def properties(self):  # {{{
    136         # ordered list of properties since vars(self) is random
    137         return ['mesh',
    138                 'mask',
    139                 'constants',
    140                 'geometry',
    141                 'initialization',
    142                 'smb',
    143                 'basalforcings',
    144                 'friction',
    145                 'rifts',
    146                 'solidearth',
    147                 'dsl',
    148                 'timestepping',
    149                 'groundingline',
    150                 'materials',
    151                 'damage',
    152                 'flowequation',
    153                 'debug',
    154                 'verbose',
    155                 'settings',
    156                 'toolkits',
    157                 'cluster',
    158                 'balancethickness',
    159                 'stressbalance',
    160                 'hydrology',
    161                 'masstransport',
    162                 'thermal',
    163                 'steadystate',
    164                 'transient',
    165                 'levelset',
    166                 'calving',
    167                 'frontalforcings',
    168                 'love',
    169                 'esa',
    170                 'sampling',
    171                 'autodiff',
    172                 'inversion',
    173                 'qmu',
    174                 'amr',
    175                 'radaroverlay',
    176                 'results',
    177                 'outputdefinition',
    178                 'miscellaneous',
    179                 'private']
    180     # }}}
     139    #}}}
    181140
    182141    def __repr__(obj):  #{{{
    … …  
    185144        #   already converted <class>.__repr__ method for examples)
    186145        #
    187 
    188         #print "Here %s the number: %d" % ("is", 37)
    189146        s = "%19s: %-22s -- %s" % ("mesh", "[%s %s]" % ("1x1", obj.mesh.__class__.__name__), "mesh properties")
    190147        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("mask", "[%s %s]" % ("1x1", obj.mask.__class__.__name__), "defines grounded and floating elements"))
    … …  
    230187        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("miscellaneous", "[%s %s]" % ("1x1", obj.miscellaneous.__class__.__name__), "miscellaneous fields"))
    231188        return s
    232     # }}}
     189    #}}}
     190
     191    def properties(self): #{{{
     192        # ordered list of properties since vars(self) is random
     193        return ['mesh',
     194                'mask',
     195                'geometry',
     196                'constants',
     197                'smb',
     198                'basalforcings',
     199                'materials',
     200                'damage',
     201                'friction',
     202                'flowequation',
     203                'timestepping',
     204                'initialization',
     205                'rifts',
     206                'dsl',
     207                'solidearth',
     208                'debug',
     209                'verbose',
     210                'settings',
     211                'toolkits',
     212                'cluster',
     213                'balancethickness',
     214                'stressbalance',
     215                'groundingline',
     216                'hydrology',
     217                'masstransport',
     218                'thermal',
     219                'steadystate',
     220                'transient',
     221                'levelset',
     222                'calving',
     223                'frontalforcings',
     224                'love',
     225                'esa',
     226                'sampling',
     227                'autodiff',
     228                'inversion',
     229                'qmu',
     230                'amr',
     231                'results',
     232                'outputdefinition',
     233                'radaroverlay',
     234                'miscellaneous',
     235                'private']
     236    #}}}
    233237
    234238    def setdefaultparameters(self, planet): #{{{
    … …  
    278282    #}}}
    279283
    280     def checkmessage(self, string):  # {{{
     284    def checkmessage(self, string): #{{{
    281285        print("model not consistent: {}".format(string))
    282286        self.private.isconsistent = False
    283287        return self
    284     # }}}
     288    #}}}
    285289    #@staticmethod
    286290
    287     def extract(self, area):  # {{{
     291    def extract(self, area): #{{{
    288292        """EXTRACT - extract a model according to an Argus contour or flag list
    289293
    … …  
    562566
    563567        return md2
    564     # }}}
    565 
    566     def extrude(md, *args):  # {{{
     568    #}}}
     569
     570    def extrude(md, *args): #{{{
    567571        """EXTRUDE - vertically extrude a 2d mesh
    568572
    … …  
    749753
    750754        return md
    751     # }}}
     755    #}}}
    752756
    753757    def collapse(md): #{{{

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

    r26023 r26301  
    66classdef sampling
    77        properties (SetAccess=public)
    8          kappa          = NaN;
    9                  tau            = 0;
    10                  beta           = NaN;
    11          phi            = 0;
    12          alpha          = 0;
    13          robin          = 0;
    14          seed           = 0;
    15                  requested_outputs      = {};
    16     end
     8                kappa             = NaN;
     9                tau               = 0;
     10                beta              = NaN;
     11                phi               = 0;
     12                alpha             = 0;
     13                robin             = 0;
     14                seed              = 0;
     15                requested_outputs = {};
     16        end
    1717        methods
    1818                function self = sampling(varargin) % {{{
    19             switch nargin
     19                        switch nargin
    2020                                case 0
    2121                                        self=setdefaultparameters(self);
    2222                                otherwise
    2323                                        error('constructor not supported');
    24             end
     24                        end
    2525                end % }}}
    26                 function list = defaultoutputs(self,md) % {{{
     26                function disp(self) % {{{
    2727
    28             list = {};
     28                        disp(sprintf('   Sampling parameters:'));
    2929
    30                 end % }}}
     30                        disp(sprintf('\n      %s','Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):'));
     31                        fielddisplay(self,'kappa','coefficient of the identity operator');
     32                        fielddisplay(self,'tau','scaling coefficient of the solution (default: 1.0)');
     33                        fielddisplay(self,'alpha','exponent in PDE operator, (default: 2.0, BiLaplacian covariance operator)');
     34
     35                        disp(sprintf('\n      %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():'));
     36                        fielddisplay(self,'robin','Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default: 0)');
     37                        fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)');
     38
     39                        disp(sprintf('\n      %s','Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):'));
     40                        fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)');
     41
     42                        disp(sprintf('\n      %s','Other parameters of stochastic sampler:'));
     43                        fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default: -1)');
     44                        fielddisplay(self,'requested_outputs','additional outputs requested (not implemented yet)');
     45
     46                end % }}}','
    3147                function self = setdefaultparameters(self) % {{{
    32            
    33             %Scaling coefficient
    34                         self.tau=1; 
    3548
    36             %Apply Robin boundary conditions
    37                         self.robin=0;   
    38            
    39             %Temporal correlation factor
    40                         self.phi=0; 
    41            
    42                         %Exponent in fraction SPDE (default=2, biLaplacian covariance
     49                        %Scaling coefficient
     50                        self.tau=1;
     51
     52                        %Apply Robin boundary conditions
     53                        self.robin=0;
     54
     55                        %Temporal correlation factor
     56                        self.phi=0;
     57
     58                        %Exponent in fraction SPDE (default: 2, biLaplacian covariance
    4359                        %operator)
    4460                        self.alpha=2; % Default
    45            
    46             %Seed for pseudorandom number generator (default -1 for random seed)
    47             self.seed=-1;
    48            
     61
     62                        %Seed for pseudorandom number generator (default: -1, for random seed)
     63                        self.seed=-1;
     64
    4965                        %default output
    5066                        self.requested_outputs={'default'};
    5167
    5268                end % }}}
    53         function md = setparameters(self,md,lc,sigma) % {{{
    54            
    55             nu = self.alpha-1;
    56             KAPPA = sqrt(8*nu)/lc;
    57             TAU = sqrt(gamma(nu)/(gamma(self.alpha)*(4*pi)*KAPPA^(2*nu)*sigma^2));
    58             md.sampling.kappa = KAPPA*ones(md.mesh.numberofvertices,1);
    59             md.sampling.tau = TAU;
    60        
     69                function list = defaultoutputs(self,md) % {{{
     70
     71                        list = {};
     72
    6173                end % }}}
    6274                function md = checkconsistency(self,md,solution,analyses) % {{{
    63            
    64             if ~ismember('SamplingAnalysis',analyses), return; end
    65                        
    66             md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
    67             md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0);
    68             md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0 1]);
    69             if(md.sampling.robin)
    70                 md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
    71             end   
    72             md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0);
    73             md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0);
    74             md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1);
    75             md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1);
    7675
    77                 end % }}}
    78                 function disp(self) % {{{
     76                        if ~ismember('SamplingAnalysis',analyses), return; end
    7977
    80             disp(sprintf('   Sampling parameters:'));
     78                        md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
     79                        md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0);
     80                        md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0 1]);
     81                        if(md.sampling.robin)
     82                                md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
     83                        end
     84                        md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0);
     85                        md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0);
     86                        md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1);
     87                        md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1);
    8188
    82                         disp(sprintf('\n      %s','Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):'));
    83             fielddisplay(self,'kappa','coefficient of the identity operator');
    84                         fielddisplay(self,'tau','scaling coefficient of the solution (default 1.0)');
    85                         fielddisplay(self,'alpha','exponent in PDE operator, (default 2.0, BiLaplacian covariance operator)');
    86          
    87                         disp(sprintf('\n      %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():'));
    88                         fielddisplay(self,'robin','Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default 0)');
    89                         fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)');         
    90            
    91             disp(sprintf('\n      %s','Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):'));
    92                         fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)');
    93            
    94                         disp(sprintf('\n      %s','Other parameters of stochastic sampler:'));
    95             fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default -1)');
    96                         fielddisplay(self,'requested_outputs','additional outputs requested (not implemented yet)');
    97  
    9889                end % }}}
    9990                function marshall(self,prefix,md,fid) % {{{
    10091
    101             WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1);
    102             WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double');
    103             WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1);
    104             WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double');
    105             WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer');
    106             WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean');
    107             WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer');
    108            
     92                        WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1);
     93                        WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double');
     94                        WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1);
     95                        WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double');
     96                        WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer');
     97                        WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean');
     98                        WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer');
     99
    109100                        %process requested outputs
    110101                        outputs = self.requested_outputs;
    … …  
    116107                        WriteData(fid,prefix,'data',outputs,'name','md.sampling.requested_outputs','format','StringArray');
    117108                end % }}}
     109                function md = setparameters(self,md,lc,sigma) % {{{
     110
     111                        nu = self.alpha-1;
     112                        KAPPA = sqrt(8*nu)/lc;
     113                        TAU = sqrt(gamma(nu)/(gamma(self.alpha)*(4*pi)*KAPPA^(2*nu)*sigma^2));
     114                        md.sampling.kappa = KAPPA*ones(md.mesh.numberofvertices,1);
     115                        md.sampling.tau = TAU;
     116
     117                end % }}}
    118118                function savemodeljs(self,fid,modelname) % {{{
    119            
    120             writejsdouble(fid,[modelname '.sampling.kappa'],self.kappa);
    121             writejsdouble(fid,[modelname '.sampling.tau'],self.tau);
    122             writejsdouble(fid,[modelname '.sampling.beta'],self.beta);
    123             writejsdouble(fid,[modelname '.sampling.phi'],self.beta);
    124             writejsdouble(fid,[modelname '.sampling.alpha'],self.alpha);
    125             writejsdouble(fid,[modelname '.sampling.robin'],self.robin);
    126             writejsdouble(fid,[modelname '.sampling.seed'],self.seed);
     119
     120                        writejsdouble(fid,[modelname '.sampling.kappa'],self.kappa);
     121                        writejsdouble(fid,[modelname '.sampling.tau'],self.tau);
     122                        writejsdouble(fid,[modelname '.sampling.beta'],self.beta);
     123                        writejsdouble(fid,[modelname '.sampling.phi'],self.beta);
     124                        writejsdouble(fid,[modelname '.sampling.alpha'],self.alpha);
     125                        writejsdouble(fid,[modelname '.sampling.robin'],self.robin);
     126                        writejsdouble(fid,[modelname '.sampling.seed'],self.seed);
    127127                        writejscellstring(fid,[modelname '.sampling.requested_outputs'],self.requested_outputs);
    128128

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

    r26023 r26301  
    11import numpy as np
     2
     3from math import *
    24
    35from checkfield import checkfield
    … …  
    1416    """
    1517
    16     def __init__(self):  # {{{
    17         self.kappa = float('NaN')
     18    def __init__(self, *args): #{{{
     19        self.kappa = np.nan
    1820        self.tau = 0
    19         self.beta = float('NaN')
    20         self.hydrostatic_adjustment = 0
     21        self.beta = np.nan
    2122        self.phi = 0
    2223        self.alpha = 0
    … …  
    2526        self.requested_outputs = []
    2627
    27         # Set defaults
    28         self.setdefaultparameters()
    29 
     28        if len(args) == 0:
     29            self.setdefaultparameters()
     30        else:
     31            raise RuntimeError('constructor not supported')
    3032    #}}}
    3133
    32     def __repr__(self):  # {{{
     34    def __repr__(self): #{{{
    3335        s = '   Sampling parameters::\n'
    34         s += '{}\n'.format(fielddisplay(self,'kappa','coefficient of the identity operator'))
    35         s += '{}\n'.format(fielddisplay(self,'tau','scaling coefficient of the solution (default 1.0)'))
    36         s += '{}\n'.format(fielddisplay(self,'alpha','exponent in PDE operator, (default 2.0, BiLaplacian covariance operator)'))
    37         s += '{}\n'.format(disp(sprintf('\n      %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():')))
    38         s += '{}\n'.format(fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)'))
    39         s += '{}\n'.format(fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)'))
    40         s += '{}\n'.format(fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default -1)'))
     36        s += '      Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):\n'
     37        s += '{}\n'.format(fielddisplay(self, 'kappa', 'coefficient of the identity operator'))
     38        s += '{}\n'.format(fielddisplay(self, 'tau', 'scaling coefficient of the solution (default: 1.0)'))
     39        s += '{}\n'.format(fielddisplay(self, 'alpha', 'exponent in PDE operator, (default: 2.0, BiLaplacian covariance operator)'))
     40
     41        s += '      Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():\n'
     42        s += '{}\n'.format(fielddisplay(self, 'robin', 'Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default: 0)'))
     43        s += '{}\n'.format(fielddisplay(self, 'beta', 'Coefficient in Robin boundary conditions (to be defined for robin = 1)'))
     44
     45        s += '      Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):\n'
     46        s += '{}\n'.format(fielddisplay(self, 'phi', 'Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)'))
     47
     48        s += '      Other parameters of stochastic sampler:\n'
     49        s += '{}\n'.format(fielddisplay(self, 'seed', 'Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default: -1)'))
    4150        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested (not implemented yet)'))
     51
    4252        return s
    4353    #}}}
    4454
    45     def defaultoutputs(self, md):  # {{{
    46         return []
    47 
    48     #}}}
    49     def setdefaultparameters(self):  # {{{
     55    def setdefaultparameters(self): #{{{
    5056        # Scaling coefficient
    5157        self.tau = 1
     58
    5259        # Apply Robin boundary conditions
    5360        self.robin = 0
     61
    5462        # Temporal correlation factor
    5563        self.phi = 0
    56         # Exponent in fraction SPDE (default=2, biLaplacian covariance operator)
     64
     65        # Exponent in fraction SPDE (default: 2, biLaplacian covariance operator)
    5766        self.alpha = 2
    58         # Seed for pseudorandom number generator (default -1 for random seed)
    59         self.alpha = -1
     67
     68        # Seed for pseudorandom number generator (default: -1, for random seed)
     69        self.seed = -1
     70
    6071        # Default output
    6172        self.requested_outputs = ['default']
     73
    6274        return self
    6375    #}}}
    6476
    65     def checkconsistency(self, md, solution, analyses):  # {{{
    66         # Early return
     77    def defaultoutputs(self, md): #{{{
     78        return []
     79    #}}}
     80
     81    def checkconsistency(self, md, solution, analyses): #{{{
    6782        if ('SamplingAnalysis' not in analyses):
    6883            return md
    6984
    70         md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>',0)
    71         md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0)
    72         md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0, 1])
     85        md = checkfield(md, 'fieldname', 'sampling.kappa', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices], '>', 0)
     86        md = checkfield(md, 'fieldname', 'sampling.tau', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     87        md = checkfield(md, 'fieldname', 'sampling.robin', 'numel', 1, 'values', [0, 1])
    7388        if md.sampling.robin:
    74             md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>',0)
    75         md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0)
    76         md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0)
    77         md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1)
    78         md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1)
     89            md = checkfield(md, 'fieldname', 'sampling.beta', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices], '>', 0)
     90        end
     91        md = checkfield(md, 'fieldname', 'sampling.phi', 'NaN', 1, 'Inf', 1, 'numel', 1, '>=', 0)
     92        md = checkfield(md, 'fieldname', 'sampling.alpha', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     93        md = checkfield(md, 'fieldname', 'sampling.seed', 'NaN', 1, 'Inf', 1, 'numel', 1)
     94        md = checkfield(md, 'fieldname', 'sampling.requested_outputs', 'stringrow', 1)
    7995
    8096        return md
    81     # }}}
     97    #}}}
    8298
    83     def marshall(self, prefix, md, fid):  # {{{
    84         WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1)
    85         WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double')
    86         WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1)
    87         WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double')
    88         WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer')
    89         WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean')
    90         WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer')
     99    def marshall(self, prefix, md, fid): #{{{
     100        WriteData(fid, prefix, 'object', self, 'fieldname', 'kappa', 'format', 'DoubleMat', 'mattype', 1)
     101        WriteData(fid, prefix, 'object', self, 'fieldname', 'tau', 'format', 'Double')
     102        WriteData(fid, prefix, 'object', self, 'fieldname', 'beta', 'format', 'DoubleMat', 'mattype', 1)
     103        WriteData(fid, prefix, 'object', self, 'fieldname', 'phi', 'format', 'Double')
     104        WriteData(fid, prefix, 'object', self, 'fieldname', 'alpha', 'format', 'Integer')
     105        WriteData(fid, prefix, 'object', self, 'fieldname', 'robin', 'format', 'Boolean')
     106        WriteData(fid, prefix, 'object', self, 'fieldname', 'seed', 'format', 'Integer')
    91107
    92     #process requested outputs
     108        # Process requested outputs
    93109        outputs = self.requested_outputs
    94110        indices = [i for i, x in enumerate(outputs) if x == 'default']
    … …  
    97113            outputs = outputscopy
    98114        WriteData(fid, prefix, 'data', outputs, 'name', 'md.sampling.requested_outputs', 'format', 'StringArray')
     115    #}}}
    99116
    100     # }}}
     117    def setparameters(self, md, lc, sigma): #{{{
     118        nu = self.alpha - 1
     119        KAPPA = pow((8 * nu), 0.5) / lc
     120        TAU = pow((math.gamma(nu) / math.gamma(self.alpha) * (4 *  np.pi) * pow(KAPPA, 2 * nu) * pow(sigma, 2)), 0.5)
     121        md.sampling.kappa = KAPPA * np.ones((md.mesh.numberofvertices, 1))
     122        md.sampling.tau = TAU
     123
     124        return md
     125    #}}}

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

    r26299 r26301  
    4545                        end
    4646                end % }}}
     47                function disp(self) % {{{
     48                        disp(sprintf('   solidearth inputs, forcings and settings:'));
     49
     50                        fielddisplay(self,'planetradius','planet radius [m]');
     51                        fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
     52                        fielddisplay(self,'requested_outputs','additional outputs requested');
     53                        fielddisplay(self,'partitionice','ice partition vector for barystatic contribution');
     54                        fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution');
     55                        fielddisplay(self,'partitionocean','ocean partition vector for barystatic contribution');
     56                        if isempty(self.external), fielddisplay(self,'external','external solution, of the type solidearthsolution'); end
     57                        self.settings.disp();
     58                        self.lovenumbers.disp();
     59                        self.rotational.disp();
     60                        if ~isempty(self.external),
     61                                self.external.disp();
     62                        end
     63
     64                end % }}}
    4765                function self = setdefaultparameters(self,planet) % {{{
    4866
    … …  
    86104                function list=defaultoutputs(self,md) % {{{
    87105                        list = {'Sealevel'};
    88                 end % }}}
    89                 function disp(self) % {{{
    90                         disp(sprintf('   solidearth inputs, forcings and settings:'));
    91 
    92                         fielddisplay(self,'planetradius','planet radius [m]');
    93                         fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
    94                         fielddisplay(self,'requested_outputs','additional outputs requested');
    95                         fielddisplay(self,'partitionice','ice partition vector for barystatic contribution');
    96                         fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution');
    97                         fielddisplay(self,'partitionocean','ocean partition vector for barystatic contribution');
    98                         if isempty(self.external), fielddisplay(self,'external','external solution, of the type solidearthsolution'); end
    99                         self.settings.disp();
    100                         self.lovenumbers.disp();
    101                         self.rotational.disp();
    102                         if ~isempty(self.external),
    103                                 self.external.disp();
    104                         end
    105 
    106106                end % }}}
    107107                function marshall(self,prefix,md,fid) % {{{

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

    r26299 r26301  
    9292
    9393                end % }}}
     94                function disp(self) % {{{
     95                        disp(sprintf('   solidearth settings:'));
     96
     97                        fielddisplay(self,'reltol','sea level change relative convergence criterion (default, NaN: not applied)');
     98                        fielddisplay(self,'abstol','sea level change absolute convergence criterion(default, NaN: not applied)');
     99                        fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
     100                        fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module (default: 1)');
     101                        fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)');
     102                        fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default: 1)');
     103                        fielddisplay(self,'isgrd','compute GRD patterns (default: 1)');
     104                        fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default: 1)');
     105                        fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)');
     106                        fielddisplay(self,'selfattraction','enables surface mass load to perturb the gravity field');
     107                        fielddisplay(self,'elastic','enables elastic deformation from surface loading');
     108                        fielddisplay(self,'viscous','enables viscous deformation from surface loading');
     109                        fielddisplay(self,'rotation','enables polar motion to feedback on the GRD fields');
     110                        fielddisplay(self,'degacc','accuracy (default: .01 deg) for numerical discretization of the Green''s functions');
     111                        fielddisplay(self,'timeacc','time accuracy (default: 1 yr) for numerical discretization of the Green''s functions');
     112                        fielddisplay(self,'grdmodel','type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)');
     113                        fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore');
     114                end % }}}
    94115                function md = checkconsistency(self,md,solution,analyses) % {{{
    95116
    … …  
    135156
    136157                end % }}}
    137                 function disp(self) % {{{
    138                         disp(sprintf('   solidearth settings:'));
    139 
    140                         fielddisplay(self,'reltol','sea level change relative convergence criterion (default, NaN: not applied)');
    141                         fielddisplay(self,'abstol','sea level change absolute convergence criterion(default, NaN: not applied)');
    142                         fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
    143                         fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module (default: 1)');
    144                         fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)');
    145                         fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default: 1)');
    146                         fielddisplay(self,'isgrd','compute GRD patterns (default: 1)');
    147                         fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default: 1)');
    148                         fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)');
    149                         fielddisplay(self,'selfattraction','enables surface mass load to perturb the gravity field');
    150                         fielddisplay(self,'elastic','enables elastic deformation from surface loading');
    151                         fielddisplay(self,'viscous','enables viscous deformation from surface loading');
    152                         fielddisplay(self,'rotation','enables polar motion to feedback on the GRD fields');
    153                         fielddisplay(self,'degacc','accuracy (default: .01 deg) for numerical discretization of the Green''s functions');
    154                         fielddisplay(self,'timeacc','time accuracy (default: 1 yr) for numerical discretization of the Green''s functions');
    155                         fielddisplay(self,'grdmodel','type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)');
    156                         fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore');
    157                 end % }}}
    158158                function marshall(self,prefix,md,fid) % {{{
    159159                        WriteData(fid,prefix,'object',self,'fieldname','reltol','name','md.solidearth.settings.reltol','format','Double');

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

    r26299 r26301  
    121121                if self.grdmodel == 2:
    122122                    raise RuntimeException('model requires a 2D mesh to run gia Ivins computations (change mesh from mesh3dsurface to mesh2d)')
    123                 else:
    124                     if self.grdmodel == 1:
    125                         raise RuntimeException('model requires a 3D surface mesh to run GRD computations (change mesh from mesh2d to mesh3dsurface)')
     123            else:
     124                if self.grdmodel == 1:
     125                    raise RuntimeException('model requires a 3D surface mesh to run GRD computations (change mesh from mesh2d to mesh3dsurface)')
    126126            if self.sealevelloading and not self.grdocean:
    127127                raise RuntimeException('solidearthsettings checkconsistency error message: need grdocean on if sealevelloading flag is set')
    … …  
    134134
    135135    def marshall(self, prefix, md, fid): #{{{
    136         WriteData(fid, prefix, 'object', self, 'fieldname', 'reltol', 'name', 'md.solidearth.settings.reltol', 'format', 'Double')
    137         WriteData(fid, prefix, 'object', self, 'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double')
    138         WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer')
    139         WriteData(fid, prefix, 'object', self, 'fieldname', 'selfattraction', 'name', 'md.solidearth.settings.selfattraction', 'format', 'Boolean')
    140         WriteData(fid, prefix, 'object', self, 'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean')
    141         WriteData(fid, prefix, 'object', self, 'fieldname', 'viscous', 'name', 'md.solidearth.settings.viscous', 'format', 'Boolean')
    142         WriteData(fid, prefix, 'object', self, 'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean')
    143         WriteData(fid, prefix, 'object', self, 'fieldname', 'grdocean', 'name', 'md.solidearth.settings.grdocean', 'format', 'Boolean')
    144         WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Integer')
    145         WriteData(fid, prefix, 'object', self, 'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer')
    146         WriteData(fid, prefix, 'object', self, 'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double')
    147         WriteData(fid, prefix, 'object', self, 'fieldname', 'timeacc', 'name', 'md.solidearth.settings.timeacc', 'format', 'Double', 'scale', md.constants.yts)
    148         WriteData(fid, prefix, 'object', self, 'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer')
    149         WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevelloading', 'name', 'md.solidearth.settings.sealevelloading', 'format', 'Integer')
    150         WriteData(fid, prefix, 'object', self, 'fieldname','isgrd', 'name', 'md.solidearth.settings.isgrd', 'format', 'Integer')
    151         WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_bp_grd', 'name', 'md.solidearth.settings.compute_bp_grd', 'format', 'Integer')
    152         WriteData(fid, prefix, 'object', self, 'fieldname', 'grdmodel', 'name', 'md.solidearth.settings.grdmodel', 'format', 'Integer')
    153         WriteData(fid, prefix, 'object', self, 'fieldname', 'cross_section_shape', 'name', 'md.solidearth.settings.cross_section_shape', 'format', 'Integer')
     136        WriteData(fid, prefix, 'object', self, 'fieldname', 'reltol', 'name', 'md.solidearth.settings.reltol', 'format', 'Double');
     137        WriteData(fid, prefix, 'object', self, 'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double');
     138        WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer');
     139        WriteData(fid, prefix, 'object', self, 'fieldname', 'selfattraction', 'name', 'md.solidearth.settings.selfattraction', 'format', 'Boolean');
     140        WriteData(fid, prefix, 'object', self, 'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean');
     141        WriteData(fid, prefix, 'object', self, 'fieldname', 'viscous', 'name', 'md.solidearth.settings.viscous', 'format', 'Boolean');
     142        WriteData(fid, prefix, 'object', self, 'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean');
     143        WriteData(fid, prefix, 'object', self, 'fieldname', 'grdocean', 'name', 'md.solidearth.settings.grdocean', 'format', 'Boolean');
     144        WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Boolean');
     145        WriteData(fid, prefix, 'object', self, 'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer');
     146        WriteData(fid, prefix, 'object', self, 'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double');
     147        WriteData(fid, prefix, 'object', self, 'fieldname', 'timeacc', 'name', 'md.solidearth.settings.timeacc', 'format', 'Double', 'scale',md.constants.yts);
     148        WriteData(fid, prefix, 'object', self, 'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer');
     149        WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevelloading', 'name', 'md.solidearth.settings.sealevelloading', 'format', 'Integer');
     150        WriteData(fid, prefix, 'object', self, 'fieldname', 'isgrd', 'name', 'md.solidearth.settings.isgrd', 'format', 'Integer');
     151        WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_bp_grd', 'name', 'md.solidearth.settings.compute_bp_grd', 'format', 'Integer');
     152        WriteData(fid, prefix, 'object', self, 'fieldname', 'grdmodel', 'name', 'md.solidearth.settings.grdmodel', 'format', 'Integer');
     153        WriteData(fid, prefix, 'object', self, 'fieldname', 'cross_section_shape', 'name', 'md.solidearth.settings.cross_section_shape', 'format', 'Integer');
    154154    #}}}
    155155

  • issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py

    r26186 r26301  
    1 def oshostname():
    2     import socket
     1def acosd(X): #{{{
     2    """ function acosd - Inverse cosine in degrees
    33
    4     return socket.gethostname()
     4    Usage:
     5        Y = acosd(X)
     6    """
     7    import numpy as np
     8
     9    return np.degrees(np.arccos(X))
     10#}}}
     11
     12def asind(X): #{{{
     13    """ function asind - Inverse sine in degrees
     14
     15    Usage:
     16        Y = asind(X)
     17    """
     18    import numpy as np
     19
     20    return np.degrees(np.arcsin(X))
     21#}}}
     22
     23def atand(X): #{{{
     24    """ function atand - Inverse tangent in degrees
     25
     26    Usage:
     27        Y = atand(X)
     28    """
     29    import numpy as np
     30
     31    return np.degrees(np.arctan(X))
     32#}}}
    533
    634
    7 def ispc():
    8     import platform
     35def atan2d(Y, X): #{{{
     36    """ function atan2d - Four-quadrant inverse tangent in degrees
    937
    10     if 'Windows' in platform.system():
    11         return True
     38    Usage:
     39        D = atan2d(Y, X)
     40    """
     41    import numpy as np
     42
     43    return np.degrees(np.arctan2(Y, X))
     44#}}}
     45
     46def det(a): #{{{
     47    if a.shape == (1, ):
     48        return a[0]
     49    elif a.shape == (1, 1):
     50        return a[0, 0]
     51    elif a.shape == (2, 2):
     52        return a[0, 0] * a[1, 1] - a[0, 1] * a[1, 0]
    1253    else:
    13         return False
     54        raise TypeError("MatlabFunc.det only implemented for shape (2, 2), not for shape %s." % str(a.shape))
     55#}}}
    1456
     57def heaviside(x): #{{{
     58    import numpy as np
    1559
    16 def ismac():
     60    y = np.zeros_like(x)
     61    y[np.nonzero(x > 0.)] = 1.
     62    y[np.nonzero(x == 0.)] = 0.5
     63
     64    return y
     65#}}}
     66
     67def ismac(): #{{{
    1768    import platform
    1869
    … …  
    2172    else:
    2273        return False
     74#}}}
    2375
    24 
    25 def strcmp(s1, s2):
    26 
    27     if s1 == s2:
    28         return True
    29     else:
    30         return False
    31 
    32 
    33 def strncmp(s1, s2, n):
    34 
    35     if s1[0:n] == s2[0:n]:
    36         return True
    37     else:
    38         return False
    39 
    40 
    41 def strcmpi(s1, s2):
    42 
    43     if s1.lower() == s2.lower():
    44         return True
    45     else:
    46         return False
    47 
    48 
    49 def strncmpi(s1, s2, n):
    50 
    51     if s1.lower()[0:n] == s2.lower()[0:n]:
    52         return True
    53     else:
    54         return False
    55 
    56 
    57 def ismember(a, s):
     76def ismember(a, s): #{{{
    5877    import numpy as np
    5978
    … …  
    6685    if not isinstance(a, np.ndarray):
    6786        b = [item in s for item in a]
    68 
    6987    else:
    7088        if not isinstance(s, np.ndarray):
    … …  
    7290            for i, item in enumerate(a.flat):
    7391                b.flat[i] = item in s
    74 
    7592        else:
    7693            b = np.in1d(a.flat, s.flat).reshape(a.shape)
    7794
    7895    return b
     96#}}}
    7997
     98def ispc(): #{{{
     99    import platform
    80100
    81 def det(a):
     101    if 'Windows' in platform.system():
     102        return True
     103    else:
     104        return False
     105#}}}
    82106
    83     if a.shape == (1, ):
    84         return a[0]
    85     elif a.shape == (1, 1):
    86         return a[0, 0]
    87     elif a.shape == (2, 2):
    88         return a[0, 0] * a[1, 1] - a[0, 1] * a[1, 0]
    89     else:
    90         raise TypeError("MatlabFunc.det only implemented for shape (2, 2), not for shape %s." % str(a.shape))
     107def oshostname(): #{{{
     108    import socket
    91109
     110    return socket.gethostname()
     111#}}}
    92112
    93 def sparse(ivec, jvec, svec, m=0, n=0, nzmax=0):
     113def sparse(ivec, jvec, svec, m=0, n=0, nzmax=0): #{{{
    94114    import numpy as np
    95115
    … …  
    105125
    106126    return a
     127#}}}
    107128
     129def strcmp(s1, s2): #{{{
     130    if s1 == s2:
     131        return True
     132    else:
     133        return False
     134#}}}
    108135
    109 def heaviside(x):
    110     import numpy as np
     136def strcmpi(s1, s2): #{{{
     137    if s1.lower() == s2.lower():
     138        return True
     139    else:
     140        return False
     141#}}}
    111142
    112     y = np.zeros_like(x)
    113     y[np.nonzero(x > 0.)] = 1.
    114     y[np.nonzero(x == 0.)] = 0.5
     143def strncmp(s1, s2, n): #{{{
     144    if s1[0:n] == s2[0:n]:
     145        return True
     146    else:
     147        return False
     148#}}}
    115149
    116     return y
     150def strncmpi(s1, s2, n): #{{{
     151    if s1.lower()[0:n] == s2.lower()[0:n]:
     152        return True
     153    else:
     154        return False
     155#}}}

  • issm/trunk-jpl/src/m/miscellaneous/PythonFuncs.py

    r25012 r26301  
    22
    33
    4 def logical_and_n(*arg):
    5 
     4def logical_and_n(*arg): #{{{
    65    if len(arg):
    76        result = arg[0]
    … …  
    98            result = np.logical_and(result, item)
    109        return result
    11 
    1210    else:
    1311        return None
     12#}}}
    1413
    15 
    16 def logical_or_n(*arg):
    17 
     14def logical_or_n(*arg): #{{{
    1815    if len(arg):
    1916        result = arg[0]
    … …  
    2118            result = np.logical_or(result, item)
    2219        return result
    23 
    2420    else:
    2521        return None
     22#}}}

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

    r25688 r26301  
    4545
    4646% Uncomment the following to make a copy of the binary input file for debugging
    47 % purposes (can be fed into scripts/ReadBin.py).
     47% purposes (can be fed into scripts/BinRead.py).
    4848% copyfile([md.miscellaneous.name '.bin'], [md.miscellaneous.name '.m.bin'])
    4949

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

    r25688 r26301  
    4646
    4747        # Uncomment the following to make a copy of the binary input file for
    48         # debugging purposes (can be fed into scripts/ReadBin.py).
     48        # debugging purposes (can be fed into scripts/BinRead.py).
    4949        # copy_cmd = "cp {}.bin {}.py.bin".format(md.miscellaneous.name, md.miscellaneous.name)
    5050        # subprocess.call(copy_cmd, shell=True)

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

    r26296 r26301  
    1515md.masstransport.spcthickness=[md.geometry.thickness;0];
    1616md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
     17
    1718%antarctica
    1819xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3;
    … …  
    2627md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
    2728posant=pos;
     29
    2830%greenland
    2931pos=find(late>60 & late<90 & longe>-75 & longe<-15);
    … …  
    9698Seustatic=md.results.TransientSolution.Sealevel;
    9799Beustatic=md.results.TransientSolution.Bed;
     100pause
    98101
    99102%eustatic + selfattraction run:

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

    r25956 r26301  
    11#Test Name: EarthSlc
    22import numpy as np
     3
     4from MatlabFuncs import *
    35
    46from gmshplanet import *
    … …  
    1214
    1315
    14 #mesh earth:
     16# Mesh earth
     17#
     18# NOTE: In MATLAB, we currently use cached mesh to account for differences in
     19# mesh generated under Linux versus under macOS
     20#
    1521md = model()
    16 md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 700.) #700 km resolution mesh
     22md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 700.) # 700 km resolution mesh
    1723
    18 #parameterize solidearth solution:
    19 #solidearth loading:
    20 md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, 1))
    21 md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, 1))
    22 md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1))
    23 md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1))
    24 md.dsl.sea_water_pressure_change_at_sea_floor = np.zeros((md.mesh.numberofvertices + 1, 1))
     24# Geometry for the bed, arbitrary thickness of 100
     25md.geometry.bed = np.zeros((md.mesh.numberofvertices, 1))
     26md.geometry.base = md.geometry.bed
     27md.geometry.thickness = 100 * np.ones((md.mesh.numberofvertices, 1))
     28md.geometry.surface = md.geometry.bed + md.geometry.thickness
    2529
    26 #antarctica
    27 late = md.mesh.lat[md.mesh.elements - 1].sum(axis=1) / 3
    28 longe = md.mesh.long[md.mesh.elements - 1].sum(axis=1) / 3
     30# Solidearth loading #{{{
     31md.masstransport.spcthickness = np.append(md.geometry.thickness, 0)
     32md.smb.mass_balance = np.zeros((md.mesh.numberofvertices, 1))
     33
     34# Antarctica
     35xe = md.mesh.x[md.mesh.elements - 1].sum(axis=1) / 3
     36ye = md.mesh.y[md.mesh.elements - 1].sum(axis=1) / 3
     37ze = md.mesh.z[md.mesh.elements - 1].sum(axis=1) / 3
     38re = pow((pow(xe, 2) + pow(ye, 2) + pow(ze, 2)), 0.5)
     39
     40late = asind(ze / re)
     41longe = atan2d(ye, xe)
    2942pos = np.where(late < -80)[0]
    30 md.solidearth.surfaceload.icethicknesschange[pos] = -100
    31 #greenland
    32 pos = np.where(np.logical_and.reduce((late > 70, late < 80, longe > -60, longe < -30)))[0]
    33 md.solidearth.surfaceload.icethicknesschange[pos] = -100
     43md.masstransport.spcthickness[md.mesh.elements[pos]] = md.masstransport.spcthickness[md.mesh.elements[pos]] - 100
     44posant = pos
    3445
    35 #elastic loading from love numbers:
    36 md.solidearth.lovenumbers = lovenumbers('maxdeg', 100)
     46# Greenland
     47pos = np.where(np.logical_and.reduce((late > 60, late < 90, longe > -75, longe < -15)))[0]
     48md.masstransport.spcthickness[md.mesh.elements[pos]] = md.masstransport.spcthickness[md.mesh.elements[pos]] - 100
     49posgre = pos
     50
     51# Elastic loading from love numbers:
     52md.solidearth.lovenumbers = lovenumbers('maxdeg', 1000)
    3753#}}}
    3854
    39 #mask: {{{
     55# Mask: {{{
    4056mask = gmtmask(md.mesh.lat, md.mesh.long)
     57oceanmask = -1 * np.ones((md.mesh.numberofvertices, 1))
     58pos = np.where(mask == 0)[0]
     59oceanmask[pos] = 1
     60
    4161icemask = np.ones((md.mesh.numberofvertices, 1))
    42 pos = np.where(mask == 0)[0]
    43 icemask[pos] = -1
    44 pos = np.where(mask[md.mesh.elements - 1].sum(axis=1) < 3)[0]
    45 icemask[md.mesh.elements[pos, :] - 1] = -1
     62icemask[md.mesh.elements[posant]] = -1
     63icemask[md.mesh.elements[posgre]] = -1
     64
    4665md.mask.ice_levelset = icemask
    47 md.mask.ocean_levelset = -icemask
     66md.mask.ocean_levelset = oceanmask
     67#}}}
    4868
    49 #make sure that the elements that have loads are fully grounded
    50 pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0]
    51 md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1
     69# Time stepping {{{
     70md.timestepping.start_time = 0
     71md.timestepping.time_step = 1
     72md.timestepping.final_time = 1
     73#}}}
    5274
    53 #make sure wherever there is an ice load, that the mask is set to ice:
    54 #pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] # TODO: Do we need to do this twice?
    55 md.mask.ice_levelset[md.mesh.elements[pos, :] - 1] = -1
    56 # }}}
     75# Masstransport
     76md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices, 1))
     77md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices, 1))
     78md.initialization.vx = np.zeros((md.mesh.numberofvertices, 1))
     79md.initialization.vy = np.zeros((md.mesh.numberofvertices, 1))
     80md.initialization.sealevel = np.zeros((md.mesh.numberofvertices, 1))
     81md.initialization.str = 0
    5782
    58 md.solidearth.settings.ocean_area_scaling = 0
     83# Materials
     84md.materials = materials('hydro')
    5985
    60 #geometry for the bed, arbitrary
    61 md.geometry.bed = -np.ones((md.mesh.numberofvertices, 1))
    62 
    63 #materials
    64 md.materials=materials('hydro')
    65 
    66 #Miscellaneous
     86# Miscellaneous
    6787md.miscellaneous.name = 'test2002'
    6888
    69 #Solution parameters
     89# Solution parameters
     90md.cluster.np = 3
    7091md.solidearth.settings.reltol = np.nan
    7192md.solidearth.settings.abstol = 1e-3
    72 md.solidearth.settings.computesealevelchange = 1
     93md.solidearth.settings.sealevelloading = 1
     94md.solidearth.settings.isgrd = 1
     95md.solidearth.settings.ocean_area_scaling = 0
     96md.solidearth.settings.grdmodel = 1
    7397
    74 #max number of iterations reverted back to 10 (i.e., the original default value)
     98# Physics
     99md.transient.issmb = 0
     100md.transient.isstressbalance = 0
     101md.transient.isthermal = 0
     102md.transient.ismasstransport = 1
     103md.transient.isslc = 1
     104md.solidearth.requested_outputs = ['Sealevel', 'Bed']
     105
     106# Max number of iterations reverted back to 10 (i.e., the original default value)
    75107md.solidearth.settings.maxiter = 10
    76108
    77 #eustatic run:
    78 md.solidearth.settings.rigid = 0
     109# Eustatic run
     110md.solidearth.settings.selfattraction = 0
    79111md.solidearth.settings.elastic = 0
    80112md.solidearth.settings.rotation = 0
    81 md = solve(md, 'Sealevelrise')
    82 Seustatic = md.results.SealevelriseSolution.Sealevel
     113md.solidearth.settings.viscous = 0
    83114
    84 #eustatic + rigid run:
    85 md.solidearth.settings.rigid = 1
     115md = solve(md, 'Transient')
     116Seustatic = md.results.TransientSolution.Sealevel
     117Beustatic = md.results.TransientSolution.Bed
     118
     119# Eustatic + selfattraction run
     120md.solidearth.settings.selfattraction = 1
    86121md.solidearth.settings.elastic = 0
    87122md.solidearth.settings.rotation = 0
    88 md = solve(md, 'Sealevelrise')
    89 Srigid = md.results.SealevelriseSolution.Sealevel
     123md.solidearth.settings.viscous = 0
     124md = solve(md, 'tr')
     125Sselfattraction = md.results.TransientSolution.Sealevel
     126Bselfattraction = md.results.TransientSolution.Bed
    90127
    91 #eustatic + rigid + elastic run:
    92 md.solidearth.settings.rigid = 1
     128# Eustatic + selfattraction + elastic run
     129md.solidearth.settings.selfattraction = 1
    93130md.solidearth.settings.elastic = 1
    94131md.solidearth.settings.rotation = 0
    95 md = solve(md, 'Sealevelrise')
    96 Selastic = md.results.SealevelriseSolution.Sealevel
     132md.solidearth.settings.viscous = 0
     133md = solve(md, 'tr')
     134Selastic = md.results.TransientSolution.Sealevel
     135Belastic = md.results.TransientSolution.Bed
    97136
    98 #eustatic + rigid + elastic + rotation run:
    99 md.solidearth.settings.rigid = 1
     137# Eustatic + selfattraction + elastic + rotation run
     138md.solidearth.settings.selfattraction = 1
    100139md.solidearth.settings.elastic = 1
    101140md.solidearth.settings.rotation = 1
    102 md = solve(md, 'Sealevelrise')
    103 Srotation = md.results.SealevelriseSolution.Sealevel
     141md.solidearth.settings.viscous = 0
     142md = solve(md, 'tr')
     143Srotation = md.results.TransientSolution.Sealevel
     144Brotation = md.results.TransientSolution.Bed
    104145
    105 #Fields and tolerances to track changes
    106 field_names = ['Eustatic', 'Rigid', 'Elastic', 'Rotation']
    107 field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13]
    108 field_values = [Seustatic, Srigid, Selastic, Srotation]
     146# Fields and tolerances to track changes
     147field_names = ['Seustatic', 'Sselfattraction', 'Selastic', 'Srotation', 'Beustatic', 'Bselfattraction', 'Belastic', 'Brotation']
     148field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13]
     149field_values = [Seustatic, Sselfattraction, Selastic, Srotation, Beustatic, Bselfattraction, Belastic, Brotation]
Note: See TracChangeset for help on using the changeset viewer.