Changeset 26059


Ignore:
Timestamp:
03/09/21 23:52:41 (4 years ago)
Author:
jdquinn
Message:

CHG: MATLAB -> Python API translations for recent SE changes; clean up

Location:
issm/trunk-jpl
Files:
3 added
40 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m

    r26047 r26059  
    3131%Initialize surface and basal forcings
    3232md.smb = initialize(md.smb,md);
    33 md.basalforcings   = initialize(md.basalforcings,md);
     33md.basalforcings = initialize(md.basalforcings,md);
    3434
    3535%Initialize ocean forcings and sealevel
    36 md.dsl= initialize(md.dsl,md);
     36md.dsl = initialize(md.dsl,md);
    3737
    3838%Deal with other boundary conditions

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

    r24260 r26059  
    33
    44def SetIceSheetBC(md):
    5     """
    6     SETICESHEETBC - Create the boundary conditions for stressbalance and thermal models for an IceSheet with no Ice Front
     5    """SETICESHEETBC - Create the boundary conditions for stressbalance and thermal models for an IceSheet with no Ice Front
    76
    8        Usage:
    9           md = SetIceSheetBC(md)
     7    Usage:
     8        md = SetIceSheetBC(md)
    109
    11        See also: SETICESHELFBC, SETMARINEICESHEETBC
     10    See also: SETICESHELFBC, SETMARINEICESHEETBC
    1211    """
    1312
    … …  
    3332    #No ice front -> do nothing
    3433
    35     #Create zeros basalforcings and smb
     34    #Initialize surface and basal forcings
    3635    md.smb.initialize(md)
    3736    md.basalforcings.initialize(md)
     37
     38    #Initialize ocean forcings and sealevel
     39    md.dsl.initialize(md)
    3840
    3941    #Deal with other boundary conditions

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

    r26047 r26059  
    77        properties (SetAccess=public)
    88
    9                 global_average_thermosteric_sea_level; %corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m)
    10                 sea_surface_height_above_geoid; %corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable (in m)
    11                 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!)
     9                global_average_thermosteric_sea_level; %Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).
     10                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).
     11                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!).
    1212
    1313        end
    … …  
    5151
    5252                        disp(sprintf('   dsl parameters:'));
    53                         fielddisplay(self,'global_average_thermosteric_sea_level','corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m)');
    54                         fielddisplay(self,'sea_surface_height_above_geoid','corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m)');
    55                         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)');
     53                        fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).');
     54                        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).');
     55                        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!).');
    5656
    5757                end % }}}

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

    r25962 r26059  
    1515
    1616    def __init__(self): #{{{
    17         self.global_average_thermosteric_sea_level_change = 0 #corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)
    18         self.sea_surface_height_change_above_geoid = float('NaN') #corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)
    19         self.sea_water_pressure_change_at_sea_floor = float('NaN') #corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)
    20         self.compute_fingerprints = 0; #do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid
     17        self.global_average_thermosteric_sea_level = np.nan # Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).
     18        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).
     19        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!).
    2120    #}}}
    2221
    2322    def __repr__(self): #{{{
    2423        s = '   dsl parameters:\n'
    25         s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)'))
    26         s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)'))
    27         s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)'))
    28         s += '{}\n'.format(fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid'))
     24        s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level', 'Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).'))
     25        s += '{}\n'.format(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).'))
     26        s += '{}\n'.format(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!).'))
    2927        return s
    3028    #}}}
    3129
    3230    def extrude(self, md): #{{{
    33         self.sea_surface_height_change_above_geoid = project3d(md, 'vector', self.sea_surface_height_change_above_geoid, 'type', 'node')
    34         self.sea_water_pressure_change_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_change_at_sea_floor, 'type', 'node')
     31        self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node')
     32        self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node')
    3533        return self
    3634    #}}}
    … …  
    4240    def checkconsistency(self, md, solution, analyses): #{{{
    4341        # Early return
    44         if 'SealevelriseAnalysis' not in analyses:
     42        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
    4543            return md
    46         if solution == 'TransientSolution' and md.transient.isslc == 0:
    47             return md
    48         md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level_change', 'NaN', 1, 'Inf', 1)
    49         md = checkfield(md, 'fieldname', 'dsl.sea_surface_height_change_above_geoid', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    50         md = checkfield(md, 'fieldname', 'dsl.sea_water_pressure_change_at_sea_floor', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    51         md = checkfield(md, 'fieldname', 'dsl.compute_fingerprints', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
    52         if self.compute_fingerprints:
    53             # Check if geodetic flag of slr is on
    54             if not md.slr.geodetic:
    55                 raise RuntimeError('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on')
     44
     45        md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level', 'NaN', 1, 'Inf', 1)
     46        md = checkfield(md, 'fieldname', 'dsl.sea_surface_height_above_geoid', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
     47        md = checkfield(md, 'fieldname', 'dsl.sea_water_pressure_at_sea_floor', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
     48
     49        if md.solidearth.settings.compute_bp_grd:
     50            md = checkfield(md, 'fieldname', dsl.sea_water_pressure_at_sea_floor, 'empty', 1)
     51
    5652        return md
    5753    # }}}
    … …  
    6056        yts = md.constants.yts
    6157        WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer')
    62         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'compute_fingerprints', 'format', 'Integer')
    63         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'global_average_thermosteric_sea_level_change', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', 1+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts)
    64         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_water_pressure_change_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts)
    65         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_surface_height_change_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts)
     58        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.
     59        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.
     60        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.
    6661    # }}}

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

    r26047 r26059  
    88
    99                modelid; %index into the multi-model ensemble, determine which field will be used.
    10                 global_average_thermosteric_sea_level; %corresponds to zostoga fields in CMIP5 archives. Specified as a temporally quantity (in m) for each ensemble.
    11                 sea_surface_height_above_geoid; %corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m) for each ensemble
    12                 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!)  for each ensemble
     10                global_average_thermosteric_sea_level; %Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.
     11                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) for each ensemble.
     12                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!) for each ensemble.
    1313
    1414        end
    … …  
    5353                        disp(sprintf('   dsl mme parameters:'));
    5454                        fielddisplay(self,'modelid','index into the multi-model ensemble, determine which field will be used.');
    55                         fielddisplay(self,'global_average_thermosteric_sea_level','corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.');
    56                         fielddisplay(self,'sea_surface_height_above_geoid','corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m) for each ensemble.');
    57                         fielddisplay(self,'sea_water_pressure_at_sea_floor','corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally quantity (in m) for each ensemble.');
     55                        fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.');
     56                        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) for each ensemble.');
     57                        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!) for each ensemble.');
    5858                end % }}}
    5959                function marshall(self,prefix,md,fid) % {{{

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

    r25688 r26059  
    1515
    1616    def __init__(self, *args): #{{{
    17         self.modelid                                        = 0 #index into the multi-model ensemble
    18         self.global_average_thermosteric_sea_level_change   = [] #corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) for each ensemble.
    19         self.sea_surface_height_change_above_geoid          = [] #corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble
    20         self.sea_water_pressure_change_at_sea_floor         = [] #corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr equivalent, not in Pa/yr!) for each ensemble
    21         self.compute_fingerprints                           = 0 #corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble
    22        
    23         nargin = len(args)
     17        self.modelid                               = 0 # Index into the multi-model ensemble
     18        self.global_average_thermosteric_sea_level = [] # Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.
     19        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) for each ensemble.
     20        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!) for each ensemble.
    2421
    25         if nargin == 0:
     22        nargs = len(args)
     23
     24        if nargs == 0:
    2625            self.setdefaultparameters()
    2726        else:
    … …  
    3231        s = '   dsl mme parameters:\n'
    3332        s += '{}\n'.format(fielddisplay(self, 'modelid', 'index into the multi-model ensemble, determines which field will be used.'))
    34         s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) for each ensemble.'))
    35         s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble.'))
    36         s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr) for each ensemble.'))
    37         s += '{}\n'.format(fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid'))
     33        s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level', 'Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.'))
     34        s += '{}\n'.format(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) for each ensemble.'))
     35        s += '{}\n'.format(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!) for each ensemble.'))
    3836        return s
    3937    #}}}
    … …  
    4442
    4543    def checkconsistency(self, md, solution, analyses): # {{{
    46         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslr):
     44        # Early return
     45        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
    4746            return md
    48         for i in range(len(self.global_average_thermosteric_sea_level_change)):
    49             md = checkfield(md, 'field', self.global_average_thermosteric_sea_level_change[i], 'NaN', 1, 'Inf', 1)
    50             md = checkfield(md, 'field', self.sea_surface_height_change_above_geoid[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    51             md = checkfield(md, 'field', self.sea_water_pressure_change_at_sea_floor[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    52         md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 1, '<=',len(self.global_average_thermosteric_sea_level_change))
    53         if self.compute_fingerprints:
    54             #check geodetic flag of slr is on
    55             if not md.solidearth.settings.computesealevelchange:
    56                 raise Exception('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on')
     47
     48        for i in range(len(self.global_average_thermosteric_sea_level)):
     49            md = checkfield(md, 'field', self.global_average_thermosteric_sea_level[i], 'NaN', 1, 'Inf', 1)
     50            md = checkfield(md, 'field', self.sea_surface_height_above_geoid[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
     51            md = checkfield(md, 'field', self.sea_water_pressure_at_sea_floor[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
     52        md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 1, '<=',len(self.global_average_thermosteric_sea_level))
     53
     54        if self.solidearth.settings.compute_bp_grd:
     55            md = checkfield(md, 'field', self.sea_water_pressure_at_sea_floor, 'empty', 1)
     56
    5757        return md
    5858    #}}}
    … …  
    6060    def marshall(self, prefix, md, fid): #{{{
    6161        WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 2, 'format', 'Integer')
    62         WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_fingerprints', 'format', 'Integer')
    6362        WriteData(fid, prefix, 'object', self, 'fieldname', 'modelid', 'format', 'Double')
    64         WriteData(fid, prefix, 'name', 'md.dsl.nummodels', 'data', len(self.global_average_thermosteric_sea_level_change), 'format', 'Integer')
    65         WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level_change', 'format', 'MatArray', 'timeseries', 1, 'timeserieslength', 2, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
    66         WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_change_at_sea_floor', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3)
    67         WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_change_above_geoid', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
     63        WriteData(fid, prefix, 'name', 'md.dsl.nummodels', 'data', len(self.global_average_thermosteric_sea_level), 'format', 'Integer')
     64        WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'MatArray', 'timeseries', 1, 'timeserieslength', 2, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
     65        WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3)
     66        WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_above_geoid', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
    6867    #}}}
    6968
    7069    def extrude(self, md): #{{{
    71         for i in range(len(self.global_average_thermosteric_sea_level_change)):
    72             self.sea_surface_height_change_above_geoid[i] = project3d(md, 'vector', self.self.sea_surface_height_change_above_geoid[i], 'type', 'node', 'layer', 1)
    73             self.sea_water_pressure_change_at_sea_floor[i] = project3d(md, 'vector', self.sea_water_pressure_change_at_sea_floor[i], 'type', 'node', 'layer', 1)
     70        for i in range(len(self.global_average_thermosteric_sea_level)):
     71            self.sea_surface_height_above_geoid[i] = project3d(md, 'vector', self.self.sea_surface_height_above_geoid[i], 'type', 'node', 'layer', 1)
     72            self.sea_water_pressure_at_sea_floor[i] = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor[i], 'type', 'node', 'layer', 1)
    7473
    7574        return self

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

    r26047 r26059  
    8080
    8181                        %need 'litho' material:
    82                         if ~isa(md.materials,'materials')
     82                        if ~isa(md.materials,'materials') | ~sum(strcmpi(md.materials.nature,'litho'))
    8383                                error('Need a ''litho'' material to run a Fourier Love number analysis');
    84                         else
    85                                 if ~sum(strcmpi(md.materials.nature,'litho')),
    86                                         error('Need a ''litho'' material to run a Fourier Love number analysis');
    87                                 end
    8884                        end
    8985

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

    r25300 r26059  
    55
    66class fourierlove(object):
    7     """
    8     Fourier Love Number class definition
     7    """FOURIERLOVE - Fourier Love Number class definition
    98
    10        Usage:
    11           fourierlove = fourierlove()
     9    Usage:
     10        fourierlove = fourierlove()
    1211    """
    1312    def __init__(self):  # {{{
    … …  
    2322        self.forcing_type = 0
    2423
    25     #set defaults
    2624        self.setdefaultparameters()
    2725    #}}}
    2826
    2927    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        #
    3032        string = '   Fourier Love class:'
    3133        string = "%s\n%s" % (string, fielddisplay(self, 'nfreq', 'number of frequencies sampled (default 1, elastic) [Hz]'))
    … …  
    7678
    7779    def checkconsistency(self, md, solution, analyses):  # {{{
     80        if 'LoveAnalysis' not in analyses:
     81            return md
    7882        md = checkfield(md, 'fieldname', 'love.nfreq', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    7983        md = checkfield(md, 'fieldname', 'love.frequencies', 'NaN', 1, 'Inf', 1, 'numel', [md.love.nfreq])
    … …  
    8993            raise RuntimeError("Degree 1 not supported for Volumetric Potential forcing. Use sh_min >= 2 for this kind of calculation.")
    9094
     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')
    9198        return md
    9299    # }}}

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

    r26047 r26059  
    8585                end % }}}
    8686                function marshall(self,prefix,md,fid) % {{{
    87                        
    88                         if (size(self.thickness==md.mesh.numberofvertices) | (self.thickness==md.mesh.numberofvertices+1)),
     87
     88                        if (size(self.thickness)==md.mesh.numberofvertices) | (size(self.thickness)==md.mesh.numberofvertices+1)),
    8989                                WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    90                         elseif (size(self.thickness==md.mesh.numberofelements) | (self.thickness==md.mesh.numberofelements+1)),
     90                        elseif (size(self.thickness)==md.mesh.numberofelements) | (size(self.thickness)==md.mesh.numberofelements+1)),
    9191                                WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    9292                        else

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

    r25688 r26059  
    5151
    5252    def checkconsistency(self, md, solution, analyses): #{{{
    53         if (solution == 'TransientSolution' and md.transient.isgia) or (solution == 'GiaSolution'):
    54             md = checkfield(md, 'fieldname', 'geometry.thickness', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    55         elif solution == 'SealevelriseSolution':
    56             md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    57         elif solution == 'LoveSolution':
    58             return
     53        if solution == 'LoveSolution':
     54            return md
    5955        else:
    6056            md = checkfield(md, 'fieldname', 'geometry.surface', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    … …  
    7167                    md.checkmessage('equality base = bed on grounded ice violated')
    7268                md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    73 
    7469        return md
    7570    # }}}
    7671
    7772    def marshall(self, prefix, md, fid): #{{{
     73        if (len(self.thickness) == md.mesh.numberofvertices) or (len(self.thickness) == md.mesh.numberofvertices + 1):
     74            WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
     75        elif (len(self.thickness) == md.mesh.numberofelements) or (len(self.thickness) == md.mesh.numberofelements + 1):
     76            WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
     77        else:
     78            raise RuntimeError('geometry thickness time series should be a vertex or element time series')
     79
    7880        WriteData(fid, prefix, 'object', self, 'fieldname', 'surface', 'format', 'DoubleMat', 'mattype', 1)
    79         WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    8081        WriteData(fid, prefix, 'object', self, 'fieldname', 'base', 'format', 'DoubleMat', 'mattype', 1)
    8182        WriteData(fid, prefix, 'object', self, 'fieldname', 'bed', 'format', 'DoubleMat', 'mattype', 1)

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

    r24793 r26059  
    77
    88class hydrologydc(object):
    9     """
    10     Hydrologydc class definition
     9    """HYDROLOGYDC class definition
    1110
    1211    Usage:
    13             hydrologydc = hydrologydc()
     12        hydrologydc = hydrologydc()
    1413    """
    1514
    … …  
    4948        self.eplflip_lock = 0
    5049
    51     #set defaults
    5250        self.setdefaultparameters()
    5351    # }}}
    5452
    5553    def __repr__(self):  # {{{
     54        # TODO:
     55        # - Convert all formatting to calls to <string>.format (see any
     56        #   already converted <class>.__repr__ method for examples)
     57        #
    5658        string = '   hydrology Dual Porous Continuum Equivalent parameters:'
    5759        string = ' - general parameters'

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

    r26047 r26059  
    3131                        %Type of stabilization to use 0:nothing 1:artificial_diffusivity
    3232                        self.stabilization     = 1;
    33       self.requested_outputs = {'default'};
     33                        self.requested_outputs = {'default'};
    3434                end % }}}
    3535                function md = checkconsistency(self,md,solution,analyses) % {{{

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

    r24213 r26059  
    55
    66class hydrologyshreve(object):
    7     """
    8     HYDROLOGYSHREVE class definition
     7    """HYDROLOGYSHREVE class definition
    98
    10        Usage:
    11           hydrologyshreve = hydrologyshreve()
     9    Usage:
     10        hydrologyshreve = hydrologyshreve()
    1211    """
    1312
    … …  
    1615        self.stabilization = 0
    1716        self.requested_outputs = []
    18     #set defaults
     17
    1918        self.setdefaultparameters()
    20 
    2119    #}}}
    2220
    2321    def __repr__(self):  # {{{
    24 
     22        # TODO:
     23        # - Convert all formatting to calls to <string>.format (see any
     24        #   already converted <class>.__repr__ method for examples)
     25        #
    2526        string = '   hydrologyshreve solution parameters:'
    2627        string = "%s\n%s" % (string, fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]'))
    … …  
    4849    def checkconsistency(self, md, solution, analyses):  # {{{
    4950        #Early return
    50         if 'HydrologyShreveAnalysis' not in analyses:
     51        if 'HydrologyShreveAnalysis' not in analyses or (solution == 'TransientSolution' and not md.transient.ishydrology):
    5152            return md
    5253
    … …  
    6162        WriteData(fid, prefix, 'object', self, 'fieldname', 'spcwatercolumn', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    6263        WriteData(fid, prefix, 'object', self, 'fieldname', 'stabilization', 'format', 'Double')
    63     #process requested outputs
     64        #process requested outputs
    6465        outputs = self.requested_outputs
    6566        indices = [i for i, x in enumerate(outputs) if x == 'default']

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

    r26047 r26059  
    2424                function list = defaultoutputs(self,md) % {{{
    2525                        list = {''};
    26                 end % }}}   
     26                end % }}}
    2727
    2828                function self = setdefaultparameters(self) % {{{
    … …  
    4646                        WriteData(fid,prefix,'name','md.hydrology.model','data',2,'format','Integer');
    4747                        WriteData(fid,prefix,'object',self,'fieldname','spcwatercolumn','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    48       outputs = self.requested_outputs;
    49       pos  = find(ismember(outputs,'default'));
    50       if ~isempty(pos),
    51         outputs(pos) = [];  %remove 'default' from outputs
    52         outputs      = [outputs defaultoutputs(self,md)]; %add defaults
    53       end
    54       WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
     48                        outputs = self.requested_outputs;
     49                        pos  = find(ismember(outputs,'default'));
     50                        if ~isempty(pos),
     51                                outputs(pos) = [];  %remove 'default' from outputs
     52                                outputs      = [outputs defaultoutputs(self,md)]; %add defaults
     53                        end
     54                        WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
    5555                end % }}}
    5656                function savemodeljs(self,fid,modelname) % {{{
    57                
     57
    5858                        writejs1Darray(fid,[modelname '.hydrology.spcwatercolumn'],self.spcwatercolumn);
    5959

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

    r26047 r26059  
    2222                sealevel            = NaN;
    2323                bottompressure      = NaN;
    24         sample              = NaN;
     24                sample              = NaN;
    2525        end
    2626        methods
    … …  
    5252                end % }}}
    5353                function self = setdefaultparameters(self) % {{{
    54 
    5554                end % }}}
    5655                function md = checkconsistency(self,md,solution,analyses) % {{{
    … …  
    9897                        if ismember('HydrologyShreveAnalysis',analyses),
    9998                                if isa(md.hydrology,'hydrologyshreve'),
    100                                         if strcmp(solution,'TransientSolution') & md.transient.ishydrology,
    101                                                 md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    102                                         end
    103                                         if strcmp(solution,'HydrologySolution'),
     99                                        if (strcmp(solution,'TransientSolution') & md.transient.ishydrology) | strcmp(solution,'HydrologySolution'),
    104100                                                md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    105101                                        end
    … …  
    111107                                end
    112108                        end
    113                         if ismember('SealevelchangeAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.isslc == 0),
     109                        if ismember('SealevelchangeAnalysis',analyses),
    114110                                if strcmp(solution,'TransientSolution') & md.transient.isslc,
    115111                                        md = checkfield(md,'fieldname','initialization.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    … …  
    134130                                                md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    135131                                        end
    136                 end
    137             end
    138             if ismember('SamplingAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.issampling == 0),
    139                 if ~isnan(md.initialization.sample)
    140                     md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    141                 end
     132                                end
     133                        end
     134                        if ismember('SamplingAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.issampling == 0),
     135                                if ~isnan(md.initialization.sample)
     136                                        md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
     137                                end
    142138                        end
    143139                end % }}}
    … …  
    159155                        fielddisplay(self,'hydraulic_potential','Hydraulic potential (for GlaDS) [Pa]');
    160156                        fielddisplay(self,'channelarea','subglacial water channel area (for GlaDS) [m2]');
    161             fielddisplay(self,'sample','Realization of a Gaussian random field');
     157                        fielddisplay(self,'sample','Realization of a Gaussian random field');
    162158                end % }}}
    163159                function marshall(self,prefix,md,fid) % {{{
    … …  
    179175                        WriteData(fid,prefix,'object',self,'fieldname','channelarea','format','DoubleMat','mattype',1);
    180176                        WriteData(fid,prefix,'object',self,'fieldname','hydraulic_potential','format','DoubleMat','mattype',1);
    181             WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1);
    182            
     177                        WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1);
     178
    183179                        if md.thermal.isenthalpy,
    184180                                if numel(self.enthalpy) <= 1,
    … …  
    195191                end % }}}
    196192                function savemodeljs(self,fid,modelname) % {{{
    197                
     193
    198194                        writejs1Darray(fid,[modelname '.initialization.vx'],self.vx);
    199195                        writejs1Darray(fid,[modelname '.initialization.vy'],self.vy);
    … …  
    210206                        writejs1Darray(fid,[modelname '.initialization.hydraulic_potential'],self.hydraulic_potential);
    211207                        writejs1Darray(fid,[modelname '.initialization.channel'],self.channelarea);
    212             writejs1Darray(fid,[modelname '.initialization.sample'],self.sample);
    213            
     208                        writejs1Darray(fid,[modelname '.initialization.sample'],self.sample);
     209
    214210                end % }}}
    215211        end

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

    r26014 r26059  
    2929        self.hydraulic_potential = np.nan
    3030        self.channelarea = np.nan
     31        self.sealevel = np.nan
     32        self.bottompressure = np.nan
    3133        self.sample = np.nan
    3234
    … …  
    6769        self.epl_head = project3d(md, 'vector', self.epl_head, 'type', 'node', 'layer', 1)
    6870        self.epl_thickness = project3d(md, 'vector', self.epl_thickness, 'type', 'node', 'layer', 1)
     71        self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node', 'layer', 1)
     72        self.bottompressure = project3d(md, 'vector', self.bottompressure, 'type', 'node', 'layer', 1)
    6973
    7074        #Lithostatic pressure by default
    … …  
    9094            md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    9195            md = checkfield(md, 'fieldname', 'initialization.vy', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
     96        if 'OceantransportAnalysis' in analyses:
     97            if solution == 'TransientSolution' and md.transient.isslc and md.transient.isoceantransport:
     98                md = checkfield(md, 'fieldname', 'initialization.bottompressure', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    9299        if 'BalancethicknessAnalysis' in analyses and solution == 'BalancethicknessSolution':
    93100            md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    … …  
    112119        if 'HydrologyShreveAnalysis' in analyses:
    113120            if hasattr(md.hydrology, 'hydrologyshreve'):
     121                if (solution == 'TransientSolution' and md.transient.ishydrology) or solution == 'HydrologySolution':
     122                    md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
     123        if 'HydrologyTwsAnalysis' in analyses:
     124            if hasattr(md.hydrology, 'hydrologytws'):
    114125                md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    115         if 'HydrologyDCInefficientAnalysis' in analyses:
    116             if hasattr(md.hydrology, 'hydrologydc'):
    117                 md = checkfield(md, 'fieldname', 'initialization.sediment_head', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    118         if 'HydrologyDCEfficientAnalysis' in analyses:
    119             if hasattr(md.hydrology, 'hydrologydc'):
    120                 if md.hydrology.isefficientlayer == 1:
    121                     md = checkfield(md, 'fieldname', 'initialization.epl_head', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    122                     md = checkfield(md, 'fieldname', 'initialization.epl_thickness', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
     126        if 'SealevelchangeAnalysis' in analyses:
     127            if solution == 'TransientSolution' and md.transient.isslc:
     128                md = checkfield(md, 'fieldname', 'initialization.sealevel', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    123129        if 'HydrologyGlaDSAnalysis' in analyses:
    124130            if hasattr(md.hydrology, 'hydrologyglads'):
    … …  
    138144        WriteData(fid, prefix, 'object', self, 'fieldname', 'vz', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
    139145        WriteData(fid, prefix, 'object', self, 'fieldname', 'pressure', 'format', 'DoubleMat', 'mattype', 1)
     146        WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1)
     147        WriteData(fid, prefix, 'object', self, 'fieldname', 'bottompressure', 'format', 'DoubleMat', 'mattype', 1)
    140148        WriteData(fid, prefix, 'object', self, 'fieldname', 'temperature', 'format', 'DoubleMat', 'mattype', 1)
    141149        WriteData(fid, prefix, 'object', self, 'fieldname', 'waterfraction', 'format', 'DoubleMat', 'mattype', 1)

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

    r26047 r26059  
    66classdef matdamageice
    77        properties (SetAccess=public)
    8                 rho_ice                    = 0.;
    9                 rho_water                  = 0.;
    10                 rho_freshwater             = 0.;
    11                 mu_water                   = 0.;
    12                 heatcapacity               = 0.;
    13                 latentheat                 = 0.;
    14                 thermalconductivity        = 0.;
    15                 temperateiceconductivity   = 0.;
     8                rho_ice                         = 0.;
     9                rho_water                       = 0.;
     10                rho_freshwater                  = 0.;
     11                mu_water                        = 0.;
     12                heatcapacity                    = 0.;
     13                latentheat                      = 0.;
     14                thermalconductivity             = 0.;
     15                temperateiceconductivity        = 0.;
    1616                effectiveconductivity_averaging = 0.;
    17                 meltingpoint               = 0.;
    18                 beta                       = 0.;
    19                 mixed_layer_capacity       = 0.;
    20                 thermal_exchange_velocity  = 0.;
    21                 rheology_B   = NaN;
    22                 rheology_n   = NaN;
    23                 rheology_law = '';
     17                meltingpoint                    = 0.;
     18                beta                            = 0.;
     19                mixed_layer_capacity            = 0.;
     20                thermal_exchange_velocity       = 0.;
     21                rheology_B                      = NaN;
     22                rheology_n                      = NaN;
     23                rheology_law                    = '';
    2424
    2525                %slc
    26                 earth_density              = 0;
     26                earth_density                   = 0;
    2727
    2828        end

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

    r25169 r26059  
    3131        self.rheology_law = ''
    3232
    33     #giaivins:
    34         self.lithosphere_shear_modulus = 0.
    35         self.lithosphere_density = 0.
    36         self.mantle_shear_modulus = 0.
    37         self.mantle_density = 0.
     33        #SLC
     34        self.earth_density = 5512  # average density of the Earth, (kg / m^3)
    3835
    39     #SLR
    40         self.earth_density = 5512  # average density of the Earth, (kg / m^3)
    4136        self.setdefaultparameters()
    4237    #}}}
    4338
    4439    def __repr__(self):  # {{{
     40        # TODO:
     41        # - Convert all formatting to calls to <string>.format (see any
     42        #   already converted <class>.__repr__ method for examples)
     43        #
    4544        string = "   Materials:"
    4645        string = "%s\n%s" % (string, fielddisplay(self, "rho_ice", "ice density [kg / m^3]"))
    … …  
    6059        string = "%s\n%s" % (string, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
    6160        string = "%s\n%s" % (string, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'"))
    62         string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
    63         string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_density", "Lithosphere density [g / cm^ - 3]"))
    64         string = "%s\n%s" % (string, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
    65         string = "%s\n%s" % (string, fielddisplay(self, "mantle_density", "Mantle density [g / cm^ - 3]"))
    6661        string = "%s\n%s" % (string, fielddisplay(self, "earth_density", "Mantle density [kg / m^ - 3]"))
    6762        return string
    … …  
    105100        self.rheology_law = 'Paterson'
    106101
    107         # GIA:
    108         self.lithosphere_shear_modulus = 6.7e10  # (Pa)
    109         self.lithosphere_density = 3.32  # (g / cm^ - 3)
    110         self.mantle_shear_modulus = 1.45e11  # (Pa)
    111         self.mantle_density = 3.34  # (g / cm^ - 3)
    112 
    113         #SLR
     102        #SLC
    114103        self.earth_density = 5512  #average density of the Earth, (kg / m^3)
    115104        return self
    … …  
    130119        md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
    131120        md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
    132 
    133         if 'GiaAnalysis' in analyses:
    134             md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1)
    135             md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1)
    136             md = checkfield(md,'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1)
    137             md = checkfield(md,'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1)
    138121
    139122        if 'SealevelriseAnalysis' in analyses:
    … …  
    161144        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2)
    162145        WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
    163         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
    164         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10.**3.)
    165         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double')
    166         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10.**3.)
    167146        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
    168147

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

    r26047 r26059  
    66classdef matenhancedice
    77        properties (SetAccess=public)
    8                 rho_ice                    = 0.;
    9                 rho_water                  = 0.;
    10                 rho_freshwater             = 0.;
    11                 mu_water                   = 0.;
    12                 heatcapacity               = 0.;
    13                 latentheat                 = 0.;
    14                 thermalconductivity        = 0.;
    15                 temperateiceconductivity   = 0.;
     8                rho_ice                         = 0.;
     9                rho_water                       = 0.;
     10                rho_freshwater                  = 0.;
     11                mu_water                        = 0.;
     12                heatcapacity                    = 0.;
     13                latentheat                      = 0.;
     14                thermalconductivity             = 0.;
     15                temperateiceconductivity        = 0.;
    1616                effectiveconductivity_averaging = 0.;
    17                 meltingpoint               = 0.;
    18                 beta                       = 0.;
    19                 mixed_layer_capacity       = 0.;
    20                 thermal_exchange_velocity  = 0.;
    21                 rheology_E   = NaN;
    22                 rheology_B   = NaN;
    23                 rheology_n   = NaN;
    24                 rheology_law = '';
     17                meltingpoint                    = 0.;
     18                beta                            = 0.;
     19                mixed_layer_capacity            = 0.;
     20                thermal_exchange_velocity       = 0.;
     21                rheology_E                      = NaN;
     22                rheology_B                      = NaN;
     23                rheology_n                      = NaN;
     24                rheology_law                    = '';
    2525
    26                 %slr
    27                 earth_density              = 0;
     26                %SLC
     27                earth_density                   = 0;
    2828
    2929        end

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

    r25171 r26059  
     1from checkfield import checkfield
    12from fielddisplay import fielddisplay
    23from project3d import project3d
    3 from checkfield import checkfield
    44from WriteData import WriteData
    55
    66
    77class matenhancedice(object):
    8     """
    9     MATICE class definition
     8    """MATICE class definition
    109
    11         Usage:
    12             matenhancedice = matenhancedice()
     10    Usage:
     11        matenhancedice = matenhancedice()
    1312    """
    1413
    … …  
    3231        self.rheology_law = ''
    3332
    34         #GIA
    35         self.lithosphere_shear_modulus = 0.
    36         self.lithosphere_density = 0.
    37         self.mantle_shear_modulus = 0.
    38         self.mantle_density = 0.
     33        #SLC
     34        self.earth_density = 0  # average density of the Earth, (kg/m^3)
    3935
    40         #SLR
    41         self.earth_density = 0  # average density of the Earth, (kg/m^3)
    4236        self.setdefaultparameters()
    4337    #}}}
    4438
    4539    def __repr__(self): #{{{
     40        # TODO:
     41        # - Convert all formatting to calls to <string>.format (see any
     42        #   already converted <class>.__repr__ method for examples)
     43        #
    4644        s = "   Materials:"
    4745        s = "%s\n%s" % (s, fielddisplay(self, "rho_ice", "ice density [kg/m^3]"))
    … …  
    6260        s = "%s\n%s" % (s, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
    6361        s = "%s\n%s" % (s, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'"))
    64         s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
    65         s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_density", "Lithosphere density [g/cm^-3]"))
    66         s = "%s\n%s" % (s, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
    67         s = "%s\n%s" % (s, fielddisplay(self, "mantle_density", "Mantle density [g/cm^-3]"))
    6862        s = "%s\n%s" % (s, fielddisplay(self, "earth_density", "Mantle density [kg/m^-3]"))
    6963
    … …  
    109103        self.rheology_law = 'Paterson'
    110104
    111     # GIA:
     105        #GIA
    112106        self.lithosphere_shear_modulus = 6.7 * 10**10  # (Pa)
    113107        self.lithosphere_density = 3.32  # (g / cm^ - 3)
    … …  
    115109        self.mantle_density = 3.34  # (g / cm^ - 3)
    116110
    117     #SLR
     111        #SLC
    118112        self.earth_density = 5512  #average density of the Earth, (kg / m^3)
    119113
    … …  
    132126        md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
    133127
    134         if 'GiaAnalysis' in analyses:
    135             md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1)
    136             md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1)
    137             md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1)
    138             md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1)
    139            
    140128        if 'SealevelriseAnalysis' in analyses:
    141129            md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', 1)
    … …  
    162150        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2)
    163151        WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
    164         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
    165         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10**3)
    166         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double')
    167         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10**3)
    168152        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
    169153    # }}}

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

    r26047 r26059  
    122122                                        self.rheology_B   = 1*1e8;
    123123                                        self.rheology_n   = 3;
    124 
    125124
    126125                                case 'litho'

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

    r25167 r26059  
    88
    99class materials(object):
    10     '''
    11     MATERIALS class definition
    12 
    13        Usage:
    14           materials = materials()
    15     '''
     10    """MATERIALS class definition
     11
     12    Usage:
     13        materials = materials()
     14    """
    1615
    1716    def __init__(self, *args): #{{{
    … …  
    3837                setattr(self, 'thermalconductivity', 0)
    3938                setattr(self, 'temperateiceconductivity', 0)
     39                setattr(self, 'effectiveconductivity_averaging', 0)
    4040                setattr(self, 'meltingpoint', 0)
    4141                setattr(self, 'beta', 0)
    … …  
    6060                setattr(self, 'rho_water', 0)
    6161                setattr(self, 'rho_freshwater', 0)
    62                 setattr(self, 'earth_density', 0)
    6362            else:
    6463                raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
     64        setattr(self, 'earth_density', 0)
    6565
    6666        #set default parameters:
    … …  
    7474                #ice density (kg/m^3)
    7575                self.rho_ice = 917.
     76
    7677                #ocean water density (kg/m^3)
    7778                self.rho_water = 1023.
     79
    7880                #fresh water density (kg/m^3)
    7981                self.rho_freshwater = 1000.
     82
    8083                #water viscosity (N.s/m^2)
    8184                self.mu_water = 0.001787
     85
    8286                #ice heat capacity cp (J/kg/K)
    8387                self.heatcapacity = 2093.
     88
    8489                #ice latent heat of fusion L (J / kg)
    8590                self.latentheat = 3.34e5
     91
    8692                #ice thermal conductivity (W/m/K)
    8793                self.thermalconductivity = 2.4
     94
    8895                #wet ice thermal conductivity (W/m/K)
    8996                self.temperateiceconductivity = 0.24
     97
     98                #computation of effective conductivity
     99                self.effectiveconductivity_averaging = 1
     100
    90101                #the melting point of ice at 1 atmosphere of pressure in K
    91102                self.meltingpoint = 273.15
     103
    92104                #rate of change of melting point with pressure (K/Pa)
    93105                self.beta = 9.8e-8
     106
    94107                #mixed layer (ice-water interface) heat capacity (J/kg/K)
    95108                self.mixed_layer_capacity = 3974.
     109
    96110                #thermal exchange velocity (ice-water interface) (m/s)
    97111                self.thermal_exchange_velocity = 1.00e-4
     112
    98113                #Rheology law: what is the temperature dependence of B with T
    99114                #available: none, paterson and arrhenius
    100115                self.rheology_law = 'Paterson'
    101             elif nat == 'litho':
    102                 #we default to a configuration that enables running GIA solutions using giacaron and / or giaivins.
     116
     117                #Rheology fields default
     118                self.rheology_B = 1e8
     119                self.rheology_n = 3
     120            elif nat == 'litho':
     121                #we default to a configuration that enables running GIA solutions using giacaron and/or giaivins.
    103122                self.numlayers = 2
     123
    104124                #center of the earth (approximation, must not be 0), then the lab (lithosphere / asthenosphere boundary) then the surface
    105125                #(with 1d3 to avoid numerical singularities)
    106126                self.radius = [1e3, 6278e3, 6378e3]
     127
    107128                self.viscosity = [1e21, 1e40]  #mantle and lithosphere viscosity (respectively) [Pa.s]
    108129                self.lame_mu = [1.45e11, 6.7e10]  # (Pa)  #lithosphere and mantle shear modulus (respectively) [Pa]
    … …  
    116137                #ice density (kg/m^3)
    117138                self.rho_ice = 917.
     139
    118140                #ocean water density (kg/m^3)
    119141                self.rho_water = 1023.
    120                 #average density of the Earth (kg/m^3)
    121                 self.earth_density = 5512
     142
    122143                #fresh water density (kg/m^3)
    123144                self.rho_freshwater = 1000.
    124145            else:
    125146                raise RuntimeError("materials setdefaultparameters error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
     147
     148            #average density of the Earth (kg/m^3)
     149            self.earth_density = 5512
    126150    #}}}
    127151
    … …  
    223247        WriteData(fid, prefix, 'name', 'md.materials.nature', 'data', naturetointeger(self.nature), 'format', 'IntMat', 'mattype', 3)
    224248        WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER, this can evolve if you have classes
    225 
    226249        for i in range(len(self.nature)):
    227250            nat = self.nature[i]
    … …  
    236259                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'thermalconductivity', 'format', 'Double')
    237260                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'temperateiceconductivity', 'format', 'Double')
     261                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'effectiveconductivity_averaging', 'format', 'Integer')
    238262                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'meltingpoint', 'format', 'Double')
    239263                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'beta', 'format', 'Double')
    … …  
    254278                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_viscosity', 'format', 'DoubleMat', 'mattype', 3)
    255279                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_mu', 'format', 'DoubleMat', 'mattype', 3)
     280                # Compute earth density compatible with our layer density distribution
     281                earth_density = 0
     282                for i in range(len(self.numlayers)):
     283                    earth_density = earth_density + (self.radius[i + 1] ** 3 - self.radius[i] ** 3) * self.density[i]
     284                earth_density = earth_density / self.radius[self.numlayers + 1] ** 3
     285                self.earth_density = earth_density
    256286            elif nat == 'hydro':
    257287                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double')
    258288                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_water', 'format', 'Double')
    259                 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
    260289                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_freshwater', 'format', 'Double')
    261290            else:
    262291                raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
     292        WriteData(fid, prefix, 'data', self.earth_density, 'name', 'md.materials.earth_density', 'format', 'Double')
    263293    #}}}
    264294

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

    r26047 r26059  
    66classdef matestar
    77        properties (SetAccess=public)
    8                 rho_ice                    = 0.;
    9                 rho_water                  = 0.;
    10                 rho_freshwater             = 0.;
    11                 mu_water                   = 0.;
    12                 heatcapacity               = 0.;
    13                 latentheat                 = 0.;
    14                 thermalconductivity        = 0.;
    15                 temperateiceconductivity   = 0.;
     8                rho_ice                         = 0.;
     9                rho_water                       = 0.;
     10                rho_freshwater                  = 0.;
     11                mu_water                        = 0.;
     12                heatcapacity                    = 0.;
     13                latentheat                      = 0.;
     14                thermalconductivity             = 0.;
     15                temperateiceconductivity        = 0.;
    1616                effectiveconductivity_averaging = 0;
    17                 meltingpoint               = 0.;
    18                 beta                       = 0.;
    19                 mixed_layer_capacity       = 0.;
    20                 thermal_exchange_velocity  = 0.;
    21                 rheology_B    = NaN;
    22                 rheology_Ec   = NaN;
    23                 rheology_Es   = NaN;
    24                 rheology_law = '';
     17                meltingpoint                    = 0.;
     18                beta                            = 0.;
     19                mixed_layer_capacity            = 0.;
     20                thermal_exchange_velocity       = 0.;
     21                rheology_B                      = NaN;
     22                rheology_Ec                     = NaN;
     23                rheology_Es                     = NaN;
     24                rheology_law                    = '';
    2525
    2626                %slc
    27                 earth_density              = 0;
     27                earth_density                   = 0;
    2828
    2929        end

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

    r25171 r26059  
    3434        self.rheology_law = ''
    3535
    36         #GIA
    37         self.lithosphere_shear_modulus = 0.
    38         self.lithosphere_density = 0.
    39         self.mantle_shear_modulus = 0.
    40         self.mantle_density = 0.
    41 
    42         #SLR
     36        #slc
    4337        self.earth_density = 0
    4438
    … …  
    6660        s = "%s\n%s" % (s, fielddisplay(self, 'rheology_Es', 'shear enhancement factor'))
    6761        s = "%s\n%s" % (s, fielddisplay(self, 'rheology_law', ['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', ''Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']))
    68         s = "%s\n%s" % (s, fielddisplay(self, 'lithosphere_shear_modulus', 'Lithosphere shear modulus [Pa]'))
    69         s = "%s\n%s" % (s, fielddisplay(self, 'lithosphere_density', 'Lithosphere density [g/cm^-3]'))
    70         s = "%s\n%s" % (s, fielddisplay(self, 'mantle_shear_modulus', 'Mantle shear modulus [Pa]'))
    71         s = "%s\n%s" % (s, fielddisplay(self, 'mantle_density', 'Mantle density [g/cm^-3]'))
    7262        s = "%s\n%s" % (s, fielddisplay(self, 'earth_density', 'Mantle density [kg/m^-3]'))
    7363
    … …  
    113103        #available: none, paterson and arrhenius
    114104        self.rheology_law = 'Paterson'
    115         # GIA
    116         self.lithosphere_shear_modulus = 6.7e10  # (Pa)
    117         self.lithosphere_density = 3.32  # (g / cm^ - 3)
    118         self.mantle_shear_modulus = 1.45e11  # (Pa)
    119         self.mantle_density = 3.34  # (g / cm^ - 3)
    120         #SLR
     105        #slc
    121106        self.earth_density = 5512  # average density of the Earth, (kg / m^3)
    122107
    … …  
    134119        md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval'])
    135120        md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
    136 
    137         if 'GiaAnalysis' in analyses:
    138             md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1)
    139             md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1)
    140             md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1)
    141             md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1)
    142121
    143122        if 'SealevelriseAnalysis' in analyses:
    … …  
    166145        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_Es', 'format', 'DoubleMat', 'mattype', 1)
    167146        WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
    168         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
    169         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 1.0e3)
    170         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double')
    171         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 1.0e3)
    172147        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
    173148    # }}}

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

    r26047 r26059  
    104104                end % }}}
    105105                function md = checkconsistency(self,md,solution,analyses) % {{{
    106                        
     106
    107107                        if strcmpi(solution,'TransientSolution') & md.transient.isslc,
    108108                                md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1);
    109                                 return;
    110                         end
    111        
    112                         md = checkfield(md,'fieldname','materials.rho_ice','>',0);
    113                         md = checkfield(md,'fieldname','materials.rho_water','>',0);
    114                         md = checkfield(md,'fieldname','materials.rho_freshwater','>',0);
    115                         md = checkfield(md,'fieldname','materials.mu_water','>',0);
    116                         md = checkfield(md,'fieldname','materials.rheology_B','>',0,'universal',1,'NaN',1,'Inf',1);
    117                         md = checkfield(md,'fieldname','materials.rheology_n','>',0,'universal',1,'NaN',1,'Inf',1);
    118                         md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval' 'NyeCO2' 'NyeH2O'});
    119                         md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]);
     109                        else
     110                                md = checkfield(md,'fieldname','materials.rho_ice','>',0);
     111                                md = checkfield(md,'fieldname','materials.rho_water','>',0);
     112                                md = checkfield(md,'fieldname','materials.rho_freshwater','>',0);
     113                                md = checkfield(md,'fieldname','materials.mu_water','>',0);
     114                                md = checkfield(md,'fieldname','materials.rheology_B','>',0,'universal',1,'NaN',1,'Inf',1);
     115                                md = checkfield(md,'fieldname','materials.rheology_n','>',0,'universal',1,'NaN',1,'Inf',1);
     116                                md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval' 'NyeCO2' 'NyeH2O'});
     117                                md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]);
    120118
    121119                end % }}}

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

    r25385 r26059  
    88
    99class matice(object):
    10     '''
    11     MATICE class definition
     10    """MATICE class definition
    1211
    1312    Usage:
    1413            matice = matice()
    15     '''
     14    """
    1615
    1716    def __init__(self): #{{{
    … …  
    3332        self.rheology_law = ''
    3433
    35         #giaivins
    36         self.lithosphere_shear_modulus = 0.
    37         self.lithosphere_density = 0.
    38         self.mantle_shear_modulus = 0.
    39         self.mantle_density = 0.
    40 
    41         #slr
     34        #slc
    4235        self.earth_density = 0
    4336        self.setdefaultparameters()
    … …  
    6255        s = "%s\n%s" % (s, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
    6356        s = "%s\n%s" % (s, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', or 'NyeH2O'"))
    64         s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
    65         s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_density", "Lithosphere density [g/cm^-3]"))
    66         s = "%s\n%s" % (s, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
    67         s = "%s\n%s" % (s, fielddisplay(self, "mantle_density", "Mantle density [g/cm^-3]"))
    6857        s = "%s\n%s" % (s, fielddisplay(self, "earth_density", "Mantle density [kg/m^-3]"))
    6958
    … …  
    10897        self.rheology_law = 'Paterson'
    10998
    110         # GIA:
    111         self.lithosphere_shear_modulus = 6.7e10 # (Pa)
    112         self.lithosphere_density = 3.32  # (g/cm^-3)
    113         self.mantle_shear_modulus = 1.45e11  # (Pa)
    114         self.mantle_density = 3.34  # (g/cm^-3)
     99        # Rheology for ice
     100        self.rheology_B = 2.1 * 1e8
     101        self.rheology_n = 3
    115102
    116103        # SLR
    … …  
    119106
    120107    def checkconsistency(self, md, solution, analyses): #{{{
    121         if solution != 'SealevelriseSolution':
     108        if solution == 'TransientSolution' and md.transient.isslc:
     109            md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
     110        else:
    122111            md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
    123112            md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0)
    … …  
    125114            md = checkfield(md, 'fieldname', 'materials.mu_water', '>', 0)
    126115            md = checkfield(md, 'fieldname', 'materials.rheology_B', '>', 0, 'universal', 1, 'NaN', 1, 'Inf', 1)
    127             md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements])
     116            md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'universal',1, 'NaN', 1, 'Inf', 1)
    128117            md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O'])
    129             md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
    130 
    131         if 'GiaAnalysis' in analyses:
    132             md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', [1])
    133             md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', [1])
    134             md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', [1])
    135             md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
    136 
    137         if 'SealevelriseAnalysis' in analyses:
    138             md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
     118            md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
    139119
    140120        return md
    … …  
    159139        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2)
    160140        WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
    161         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
    162         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10.**3.)
    163         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double')
    164         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10.**3.)
    165141        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
    166142    #}}}

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

    r26047 r26059  
    4444                love                     = 0;
    4545                esa              = 0;
    46         sampling         = 0;
     46                sampling         = 0;
    4747
    4848                autodiff         = 0;
    … …  
    156156                        if isa(md.frontalforcings,'double');
    157157                                if(isprop('meltingrate',md.calving) & ~isnan(md.calving.meltingrate))
    158                                         disp('Warning: md.calving.meltingrate is now in md.frontalforcings');
     158                gia                     disp('Warning: md.calving.meltingrate is now in md.frontalforcings');
    159159                                end
    160160                                md.frontalforcings=frontalforcings(md.calving);
    … …  
    187187                                        end
    188188                                end
    189             end
    190             %2021 February 17
     189                        end
     190                        %2021 February 17
    191191                        if isa(md.sampling,'double'); md.sampling=sampling(); end
    192192                end% }}}
    … …  
    242242                        md.love             = fourierlove();
    243243                        md.esa              = esa();
    244             md.sampling         = sampling();
     244                        md.sampling         = sampling();
    245245                        md.autodiff         = autodiff();
    246246                        md.inversion        = inversion();
    … …  
    615615                                                elseif fieldsize(1)==numberofelements1
    616616                                                        md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:);
    617                         elseif (fieldsize(1)==numberofelements1+1)
    618                             md2.(model_fields{i}).(object_fields{j})=[field(pos_elem,:); field(end,:)];
     617                                                elseif (fieldsize(1)==numberofelements1+1)
     618                                                        md2.(model_fields{i}).(object_fields{j})=[field(pos_elem,:); field(end,:)];
    619619                                                end
    620620                                        end
    … …  
    628628                                        elseif fieldsize(1)==numberofelements1
    629629                                                md2.(model_fields{i})=field(pos_elem,:);
    630                     elseif (fieldsize(1)==numberofelements1+1)
     630                                        elseif (fieldsize(1)==numberofelements1+1)
    631631                                                md2.(model_fields{i})=[field(pos_elem,:); field(end,:)];
    632632                                        end
    … …  
    775775                                                % loop over time steps
    776776                                                for p=1:length(md1.results.(solutionfields{i}))
    777                                                     current = md1.results.(solutionfields{i})(p);
    778                                                     solutionsubfields=fields(current);
    779                                                     for j=1:length(solutionsubfields),
     777                                                        current = md1.results.(solutionfields{i})(p);
     778                                                        solutionsubfields=fields(current);
     779                                                        for j=1:length(solutionsubfields),
    780780                                                        field=md1.results.(solutionfields{i})(p).(solutionsubfields{j});
    781781                                                        if length(field)==numberofvertices1,
    782                                                             md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_node);
     782                                                                md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_node);
    783783                                                        elseif length(field)==numberofelements1,
    784                                                             md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_elem);
     784                                                                md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_elem);
    785785                                                        else
    786                                                             md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field;
     786                                                                md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field;
    787787                                                        end
    788                                                     end
     788                                                        end
    789789                                                end
    790790                                        else
    … …  
    15011501                        disp(sprintf('%19s: %-22s -- %s','esa'             ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution'));
    15021502                        disp(sprintf('%19s: %-22s -- %s','love'            ,['[1x1 ' class(self.love) ']'],'parameters for love solution'));
    1503             disp(sprintf('%19s: %-22s -- %s','sampling'        ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler'));
     1503                        disp(sprintf('%19s: %-22s -- %s','sampling'        ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler'));
    15041504                        disp(sprintf('%19s: %-22s -- %s','autodiff'        ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters'));
    15051505                        disp(sprintf('%19s: %-22s -- %s','inversion'       ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods'));

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

    r26013 r26059  
    5555from steadystate import steadystate
    5656from transient import transient
    57 from giaivins import giaivins
    58 from giamme import giamme
    5957from esa import esa
    6058from autodiff import autodiff
    … …  
    114112        self.calving = None
    115113        self.frontalforcings = None
    116         self.gia = None
    117114        self.love = None
    118115        self.esa = None
    … …  
    172169                'calving',
    173170                'frontalforcings',
    174                 'gia',
    175171                'love',
    176172                'esa',
    … …  
    225221        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("calving", "[%s %s]" % ("1x1", obj.calving.__class__.__name__), "parameters for calving"))
    226222        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("frontalforcings", "[%s %s]" % ("1x1", obj.frontalforcings.__class__.__name__), "parameters for frontalforcings"))
    227         s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("gia", "[%s %s]" % ("1x1", obj.gia.__class__.__name__), "parameters for gia solution"))
    228223        s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("esa", "[%s %s]" % ("1x1", obj.esa.__class__.__name__), "parameters for elastic adjustment solution"))
    229224        s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("sampling", "[%s %s]" % ("1x1", obj.sampling.__class__.__name__), "parameters for stochastic sampler"))
    … …  
    272267        self.calving = calving()
    273268        self.frontalforcings = frontalforcings()
    274         self.gia = giamme()
    275269        self.love = fourierlove()
    276270        self.esa = esa()
    … …  
    852846            md.initialization.watercolumn = project2d(md, md.initialization.watercolumn, 1)
    853847
    854         # giaivins
    855         if md.gia.__class__.__name__ == 'giaivins':
    856             if not np.isnan(md.gia.mantle_viscosity).all():
    857                 md.gia.mantle_viscosity = project2d(md, md.gia.mantle_viscosity, 1)
    858             if not np.isnan(md.gia.lithosphere_thickness).all():
    859                 md.gia.lithosphere_thickness = project2d(md, md.gia.lithosphere_thickness, 1)
    860 
    861848        # elementstype
    862849        if not np.isnan(md.flowequation.element_equation).all():

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

    r26047 r26059  
    9191                                md= slm.icecaps{i};
    9292                                if ~isempty(md.solidearth.external),
    93                                         error(sprintf('cannot run external forcings on an ice sheet when runing a coupling earth/ice sheet model');
    94                                 end
    95 
    96                         end
    97                         %make sure that we have the rigth grd model for computing ouf sealevel patterns:
     93                                        error('cannot run external forcings on an ice sheet when running a coupling earth/ice sheet model');
     94                                end
     95
     96                        end
     97                        %make sure that we have the right grd model for computing out sealevel patterns:
    9898                        for i=1:length(slm.icecaps),
    9999                                md= slm.icecaps{i};
    100100                                if md.solidearth.settings.grdmodel~=0
    101                                         error(sprintf('sealevelmodel checkconsistenty error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap %i',i));
     101                                        error(sprintf('sealevelmodel checkconsistency error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap %i',i));
    102102                                end
    103103                        end

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

    r25688 r26059  
    111111            if np.nonzero(md.slr.steric_rate - slm.earth.slr.steric_rate[slm.earth.slr.transitions[i]]) != []:
    112112                raise Exception('steric rate on ice cap {} is not the same as for the earth'.format(md.miscellaneous.name))
     113
     114        # Make sure grd is the same everywhere
     115        for i in range(len(slm.icecaps)):
     116            md = slm.icecaps[i]
     117            if md.solidearthsettings.isgrd != slm.earth.solidearthsettings.isgrd:
     118                raise RuntimeError('isgrd on ice cap {} is not the same as for the earth\n'.format(md.miscellaneous.name))
     119        # Make sure that there is no solid earth external forcing on the basins
     120        for i in range(len(slm.icecaps)):
     121            md = slm.icecaps[i]
     122            if md.solidearth.external:
     123                raise RuntimeError('cannot run external forcings on an ice sheet when running a coupling earth/ice sheet model')
     124        # Make sure that we have the right grd model for computing our sealevel patterns
     125        for i in range(len(slm.icecaps)):
     126            md = slm.icecaps[i]
     127            if md.solidearth.settings.grdmodel != 0:
     128                raise RuntimeError('sealevelmodel checkconsistency error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap {}'.format(i))
    113129    #}}}
    114130

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

    r25960 r26059  
    2020                end % }}}
    2121                function self = setdefaultparameters(self) % {{{
    22                
     22
    2323                        icethicknesschange=[];
    2424                        waterheightchange=[];
    2525                        otherchange=[];
    26                
     26
    2727                end % }}}
    28                         function md = checkconsistency(self,md,solution,analyses) % {{{
     28                function md = checkconsistency(self,md,solution,analyses) % {{{
    2929
    3030                        if ~ismember('SealevelchangeAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.isslc==0),

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

    r26047 r26059  
    88                issmb             = 0;
    99                ismasstransport   = 0;
    10                 isoceantransport   = 0;
     10                isoceantransport  = 0;
    1111                isstressbalance   = 0;
    1212                isthermal         = 0;
    … …  
    1616                ismovingfront     = 0;
    1717                ishydrology       = 0;
    18         issampling        = 0;
     18                issampling        = 0;
    1919                isslc             = 0;
    2020                amr_frequency     = 0;
    … …  
    4444                        self.ismovingfront   =0;
    4545                        self.ishydrology     = 0;
    46             self.issampling      = 0;
     46                        self.issampling      = 0;
    4747                        self.isslc           = 0;
    4848                        self.isoceancoupling = 0;
    … …  
    6565                        self.ismovingfront   = 0;
    6666                        self.ishydrology     = 0;
    67             self.issampling      = 0;
     67                        self.issampling      = 0;
    6868                        self.isslc           = 0;
    6969                        self.isoceancoupling = 0;
    … …  
    9898                        md = checkfield(md,'fieldname','transient.isslc','numel',[1],'values',[0 1]);
    9999                        md = checkfield(md,'fieldname','transient.isoceancoupling','numel',[1],'values',[0 1]);
    100             md = checkfield(md,'fieldname','transient.issampling','numel',[1],'values',[0 1]); 
     100                        md = checkfield(md,'fieldname','transient.issampling','numel',[1],'values',[0 1]); 
    101101                        md = checkfield(md,'fieldname','transient.amr_frequency','numel',[1],'>=',0,'NaN',1,'Inf',1);
    102102
    … …  
    121121                        fielddisplay(self,'ismovingfront','indicates whether a moving front capability is used in the transient');
    122122                        fielddisplay(self,'ishydrology','indicates whether an hydrology model is used');
    123             fielddisplay(self,'issampling','indicates whether sampling is used in the transient')
     123                        fielddisplay(self,'issampling','indicates whether sampling is used in the transient')
    124124                        fielddisplay(self,'isslc','indicates whether a sea-level change solution is used in the transient');
    125125                        fielddisplay(self,'isoceancoupling','indicates whether a coupling with an ocean model is used in the transient');
    … …  
    139139                        WriteData(fid,prefix,'object',self,'fieldname','ishydrology','format','Boolean');
    140140                        WriteData(fid,prefix,'object',self,'fieldname','ismovingfront','format','Boolean');
    141             WriteData(fid,prefix,'object',self,'fieldname','issampling','format','Boolean');
     141                        WriteData(fid,prefix,'object',self,'fieldname','issampling','format','Boolean');
    142142                        WriteData(fid,prefix,'object',self,'fieldname','isslc','format','Boolean');
    143143                        WriteData(fid,prefix,'object',self,'fieldname','isoceancoupling','format','Boolean');
    … …  
    165165                        writejsdouble(fid,[modelname '.trans.ismovingfront'],self.ismovingfront);
    166166                        writejsdouble(fid,[modelname '.trans.ishydrology'],self.ishydrology);
    167             writejsdouble(fid,[modelname '.trans.issampling'],self.issampling);
     167                        writejsdouble(fid,[modelname '.trans.issampling'],self.issampling);
    168168                        writejsdouble(fid,[modelname '.trans.isslc'],self.isslc);
    169169                        writejsdouble(fid,[modelname '.trans.isoceancoupling'],self.isoceancoupling);

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

    r26011 r26059  
    1414        self.issmb = 0
    1515        self.ismasstransport = 0
     16        self.isoceantransport = 0
    1617        self.isstressbalance = 0
    1718        self.isthermal = 0
    1819        self.isgroundingline = 0
    19         self.isgia = 0
    2020        self.isesa = 0
    2121        self.isdamageevolution = 0
    … …  
    2424        self.issampling = 0
    2525        self.isslc = 0
    26         self.iscoupler = 0
    2726        self.amr_frequency = 0
    2827        self.isoceancoupling = 0
    … …  
    3635        s += '{}\n'.format(fielddisplay(self, 'issmb', 'indicates if a surface mass balance solution is used in the transient'))
    3736        s += '{}\n'.format(fielddisplay(self, 'ismasstransport', 'indicates if a masstransport solution is used in the transient'))
     37        s += '{}\n'.format(fielddisplay(self, 'isoceantransport', 'indicates whether an ocean masstransport solution is used in the transient'))
    3838        s += '{}\n'.format(fielddisplay(self, 'isstressbalance', 'indicates if a stressbalance solution is used in the transient'))
    3939        s += '{}\n'.format(fielddisplay(self, 'isthermal', 'indicates if a thermal solution is used in the transient'))
    4040        s += '{}\n'.format(fielddisplay(self, 'isgroundingline', 'indicates if a groundingline migration is used in the transient'))
    41         s += '{}\n'.format(fielddisplay(self, 'isgia', 'indicates if a postglacial rebound is used in the transient'))
    4241        s += '{}\n'.format(fielddisplay(self, 'isesa', 'indicates whether an elastic adjustment model is used in the transient'))
    4342        s += '{}\n'.format(fielddisplay(self, 'isdamageevolution', 'indicates whether damage evolution is used in the transient'))
    … …  
    4746        s += '{}\n'.format(fielddisplay(self, 'isslc', 'indicates if a sea level change solution is used in the transient'))
    4847        s += '{}\n'.format(fielddisplay(self, 'isoceancoupling', 'indicates whether coupling with an ocean model is used in the transient'))
    49         s += '{}\n'.format(fielddisplay(self, 'iscoupler', 'indicates whether different models are being run with need for coupling'))
    5048        s += '{}\n'.format(fielddisplay(self, 'amr_frequency', 'frequency at which mesh is refined in simulations with multiple time_steps'))
    5149        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'list of additional outputs requested'))
    … …  
    6361        self.issmb = 0
    6462        self.ismasstransport = 0
     63        self.isoceantransport = 0
    6564        self.isstressbalance = 0
    6665        self.isthermal = 0
    6766        self.isgroundingline = 0
    68         self.isgia = 0
    6967        self.isesa = 0
    7068        self.isdamageevolution = 0
    … …  
    7472        self.isslc = 0
    7573        self.isoceancoupling = 0
    76         self.iscoupler = 0
    7774        self.amr_frequency = 0
    7875
    … …  
    8683        self.issmb = 1
    8784        self.ismasstransport = 1
     85        self.isoceantransport = 0
    8886        self.isstressbalance = 1
    8987        self.isthermal = 1
    9088        self.isgroundingline = 0
    91         self.isgia = 0
    9289        self.isesa = 0
    9390        self.isdamageevolution = 0
    … …  
    9794        self.isslc = 0
    9895        self.isoceancoupling = 0
    99         self.iscoupler = 0
    10096        self.amr_frequency = 0
    10197
    … …  
    112108        md = checkfield(md, 'fieldname', 'transient.issmb', 'numel', [1], 'values', [0, 1])
    113109        md = checkfield(md, 'fieldname', 'transient.ismasstransport', 'numel', [1], 'values', [0, 1])
     110        md = checkfield(md, 'fieldname', 'transient.isocealtransport', 'numel', [1], 'values', [0, 1])
    114111        md = checkfield(md, 'fieldname', 'transient.isstressbalance', 'numel', [1], 'values', [0, 1])
    115112        md = checkfield(md, 'fieldname', 'transient.isthermal', 'numel', [1], 'values', [0, 1])
    116113        md = checkfield(md, 'fieldname', 'transient.isgroundingline', 'numel', [1], 'values', [0, 1])
    117         md = checkfield(md, 'fieldname', 'transient.isgia', 'numel', [1], 'values', [0, 1])
    118114        md = checkfield(md, 'fieldname', 'transient.isesa', 'numel', [1], 'values', [0, 1])
    119115        md = checkfield(md, 'fieldname', 'transient.isdamageevolution', 'numel', [1], 'values', [0, 1])
    … …  
    123119        md = checkfield(md, 'fieldname', 'transient.isslc', 'numel', [1], 'values', [0, 1])
    124120        md = checkfield(md, 'fieldname', 'transient.isoceancoupling', 'numel', [1], 'values', [0, 1])
    125         md = checkfield(md, 'fieldname', 'transient.iscoupler', 'numel', [1], 'values', [0, 1])
    126121        md = checkfield(md, 'fieldname', 'transient.amr_frequency', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1)
    127122
    … …  
    136131        WriteData(fid, prefix, 'object', self, 'fieldname', 'issmb', 'format', 'Boolean')
    137132        WriteData(fid, prefix, 'object', self, 'fieldname', 'ismasstransport', 'format', 'Boolean')
     133        WriteData(fid, prefix, 'object', self, 'fieldname', 'isoceantransport', 'format', 'Boolean')
    138134        WriteData(fid, prefix, 'object', self, 'fieldname', 'isstressbalance', 'format', 'Boolean')
    139135        WriteData(fid, prefix, 'object', self, 'fieldname', 'isthermal', 'format', 'Boolean')
    140136        WriteData(fid, prefix, 'object', self, 'fieldname', 'isgroundingline', 'format', 'Boolean')
    141         WriteData(fid, prefix, 'object', self, 'fieldname', 'isgia', 'format', 'Boolean')
    142137        WriteData(fid, prefix, 'object', self, 'fieldname', 'isesa', 'format', 'Boolean')
    143138        WriteData(fid, prefix, 'object', self, 'fieldname', 'isdamageevolution', 'format', 'Boolean')
    … …  
    147142        WriteData(fid, prefix, 'object', self, 'fieldname', 'isslc', 'format', 'Boolean')
    148143        WriteData(fid, prefix, 'object', self, 'fieldname', 'isoceancoupling', 'format', 'Boolean')
    149         WriteData(fid, prefix, 'object', self, 'fieldname', 'iscoupler', 'format', 'Boolean')
    150144        WriteData(fid, prefix, 'object', self, 'fieldname', 'amr_frequency', 'format', 'Integer')
    151145

  • issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m

    r26047 r26059  
    33%
    44%   Usage:
    5 %      ismodelselfconsistent(md),
     5%      ismodelselfconsistent(md);
    66
    77%initialize consistency as true

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

    r26021 r26059  
    11def ismodelselfconsistent(md):  #{{{
    2     '''
    3     ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
     2    """ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
    43
    5        Usage:
    6           ismodelselfconsistent(md),
    7     '''
     4    Usage:
     5        ismodelselfconsistent(md)
     6    """
    87
    98    #initialize consistency as true
    … …  
    5352    elif solutiontype == 'MasstransportSolution':
    5453        analyses = ['MasstransportAnalysis']
     54    elif solutiontype == 'OceantransportSolution':
     55        analyses = ['OceantransportAnalysis']
    5556    elif solutiontype == 'BalancethicknessSolution':
    5657        analyses = ['BalancethicknessAnalysis']
    … …  
    7273        analyses = ['EsaAnalysis']
    7374    elif solutiontype == 'TransientSolution':
    74         analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'SealevelriseAnalysis']
     75        analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'OceantransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyShreveAnalysis', 'HydrologyTwsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'SealevelriseAnalysis']
    7576    elif solutiontype == 'SealevelriseSolution':
    7677        analyses = ['SealevelriseAnalysis']
    7778    elif solutiontype == 'HydrologySolution':
    78         analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis']
     79        analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'HydrologyGladsAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyTwsAnalysis']
    7980    elif 'DamageEvolutionSolution':
    8081        analyses = ['DamageEvolutionAnalysis']

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

    r26047 r26059  
    1010%   - 'Stressbalance'      or 'sb'
    1111%   - 'Masstransport'      or 'mt'
    12 %   - 'Oceantransport' or 'oceant'
     12%   - 'Oceantransport'     or 'oceant'
    1313%   - 'Thermal'            or 'th'
    1414%   - 'Steadystate'        or 'ss'

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

    r26013 r26059  
    11from datetime import datetime
    22import os
    3 import shutil
    43
    54from ismodelselfconsistent import ismodelselfconsistent
    … …  
    2221    - 'Stressbalance'      or 'sb'
    2322    - 'Masstransport'      or 'mt'
     23    - 'Oceantransport'     or 'oceant'
    2424    - 'Thermal'            or 'th'
    2525    - 'Steadystate'        or 'ss'
    … …  
    3232    - 'DamageEvolution'    or 'da'
    3333    - 'Gia'                or 'gia'
     34    - 'Love'               or 'lv'
    3435    - 'Esa'                or 'esa'
    35     - 'Sealevelchange'     or 'slc'
    36     - 'Love'               or 'lv'
    3736    - 'Sampling'           or 'smp'
    3837
    … …  
    5554    elif solutionstring.lower() == 'mt' or solutionstring.lower() == 'masstransport':
    5655        solutionstring = 'MasstransportSolution'
     56    elif solutionstring.lower() == 'oceant' or solutionstring.lower() == 'oceantransport':
     57        solutionstring = 'OceantransportSolution'
    5758    elif solutionstring.lower() == 'th' or solutionstring.lower() == 'thermal':
    5859        solutionstring = 'ThermalSolution'
    … …  
    7980    elif solutionstring.lower() == 'esa':
    8081        solutionstring = 'EsaSolution'
    81     elif solutionstring.lower() == 'slc' or solutionstring.lower() == 'sealevelchange':
    82         solutionstring = 'SealevelchangeSolution'
    8382    elif solutionstring.lower() == 'smp' or solutionstring.lower() == 'sampling':
    8483        solutionstring = 'SamplingSolution'
    … …  
    151150    if md.settings.waitonlock > 0:
    152151        islock = waitonlock(md)
    153         if islock == 0:  # no results to be loaded
     152        if islock == 0: # no results to be loaded
    154153            print('The results must be loaded manually with md = loadresultsfromcluster(md).')
    155         else:  # load results
     154        else: # load results
    156155            if md.verbose.solution:
    157156                print('loading results from cluster')

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

    r26047 r26059  
    11function slm=solveslm(slm,solutionstringi,varargin)
    2 %SOLVESLR - apply solution sequence for this sealevel model
     2%SOLVESLM - apply solution sequence for this sealevel model
    33%
    44%   Usage:
    … …  
    4343cluster=slm.cluster;
    4444batch=0;
    45 %now, go through icecaps, glacies and earth, and upload all the data independently:
     45%now, go through icecaps, glaciers and earth, and upload all the data independently:
    4646disp('solving ice caps first');
    4747for i=1:length(slm.icecaps),
    … …  
    5151slm.earth=solve(slm.earth,solutionstringi,'batch','yes');
    5252
    53 %Firs, build a runtime name that is unique
     53%First, build a runtime name that is unique
    5454c=clock;
    5555slm.private.runtimename=sprintf('%s-%02i-%02i-%04i-%02i-%02i-%02i-%i',slm.miscellaneous.name,c(2),c(3),c(1),c(4),c(5),floor(c(6)),feature('GetPid'));

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

    r26047 r26059  
    2222md.timestepping.start_time=-2400000; %4,800 kyr :: EVALUATION TIME
    2323md.timestepping.time_step= 2400000; %2,400 kyr :: EVALUATION TIME
    24 % to get rid of default final_time: make sure final_time>start_time
    25 md.timestepping.final_time=2400000; %2,500 kyr
     24% to get rid of default final_time, make sure final_time > start_time
     25md.timestepping.final_time=2400000; %2,400 kyr
    2626md.masstransport.spcthickness=[...
    2727        [md.geometry.thickness; 0],...
    … …  
    4646md.verbose=verbose('11111111111');
    4747md.verbose.solver=0;
    48 md=solve(md,'tr');
     48md=solve(md,'Transient');
    4949
    5050%Fields and tolerances to track changes

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

    r25300 r26059  
    44import numpy as np
    55
    6 from giaivins import giaivins
     6from materials import *
    77from model import *
    88from parameterize import *
    … …  
    1717md = parameterize(md, '../Par/SquareSheetConstrained.py')
    1818
    19 #GIA
    20 md.gia = giaivins()
    21 md.gia.lithosphere_thickness = 100. * np.ones(md.mesh.numberofvertices) # in km
    22 md.gia.mantle_viscosity = 1.0e21 * np.ones(md.mesh.numberofvertices) # in Pa.s
    23 md.materials.lithosphere_shear_modulus = 6.7e10 # in Pa
    24 md.materials.lithosphere_density = 3.32 # in g/cm^3
    25 md.materials.mantle_shear_modulus = 1.45e11 # in Pa
    26 md.materials.mantle_density = 3.34 # in g/cm^3
     19# GIA Ivins, 2 layer model
     20md.solidearth.settings.grdmodel = 2
    2721
    28 #Indicate what you want to compute
    29 md.gia.cross_section_shape = 1 # for square-edged x-section
     22md.materials = materials('litho','ice')
     23md.materials.numlayers = 2;
     24md.materials.radius = [10, 6271e3, 6371e3]
     25md.materials.density = [3.34e3, 3.32e3]
     26md.materials.lame_mu = [1.45e11, 6.7e10]
     27md.materials.viscosity = [1e21, 0]
     28md.initialization.sealevel = np.zeros(md.mesh.numberofvertices)
     29md.solidearth.settings.cross_section_shape = 1 # for square-edged x-section
     30md.solidearth.settings.computesealevelchange = 0 # do not compute sea level, only deformation
     31md.solidearth.requested_outputs = ['Sealevel', 'SealevelUGrd']
    3032
    31 #Define loading history (see test2001.m for the description)
    32 md.timestepping.start_time = 2400000  # 2, 400 kyr
    33 md.timestepping.final_time = 2500000  # 2, 500 kyr
    34 md.geometry.thickness = np.array([
    35     np.append(md.geometry.thickness * 0.0, 0.0),
    36     np.append(md.geometry.thickness / 2.0, 0.1),
    37     np.append(md.geometry.thickness, 0.2),
    38     np.append(md.geometry.thickness, 1.0),
    39     np.append(md.geometry.thickness, md.timestepping.start_time)
     33# Loading history
     34md.timestepping.start_time = -2400000  # 4,800 kyr :: EVALUATION TIME
     35md.timestepping.time_step = 2400000  # 2,400 kyr :: EVALUATION TIME
     36# To get rid of default final_time, make sure final_time > start_time
     37md.timestepping.final_time = 2400000  # 2,400 kyr
     38md.masstransport.spcthickness np.array([
     39    np.append(md.geometry.thickness, 0),
     40    np.append(md.geometry.thickness, 2400000)
    4041    ]).T
    4142
     43# Geometry at 0 initially
     44md.geometry.thickness = np.zeros(md.mesh.numberofvertices)
     45md.geometry.surface = np.zeros(md.mesh.numberofvertices)
     46md.geometry.base = np.zeros(md.mesh.numberofvertices)
     47
     48# Physics
     49md.transient.issmb = 0
     50md.transient.isstressbalance = 0
     51md.transient.isthermal = 0
     52md.transient.ismasstransport = 0
     53md.transient.isslc = 0
     54
    4255#Solve for GIA deflection
     56md.cluster = generic('name', gethostname(), 'np', 1)
    4357md.cluster = generic('name', gethostname(), 'np', 3)
    44 md.verbose = verbose('1111111')
    45 md = solve(md, 'Gia')
     58md.verbose = verbose('11111111111')
     59md.verbose.solver = 0
     60md = solve(md, 'Transient')
    4661
    4762#Fields and tolerances to track changes
    48 field_names = ['UGia', 'UGiaRate']
    49 field_tolerances = [1e-13, 1e-13]
    50 field_values = [md.results.GiaSolution.UGia, md.results.GiaSolution.UGiaRate]
     63field_names = ['UGrd']
     64field_tolerances = [1e-13]
     65field_values = [md.results.TransientSolution[0].SealevelUGrd]
Note: See TracChangeset for help on using the changeset viewer.