Changeset 26303


Ignore:
Timestamp:
06/08/21 00:40:42 (4 years ago)
Author:
jdquinn
Message:

CHG: More translation; test2003.py

Location:
issm/trunk-jpl
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/m/classes/timestepping.m

    r26208 r26303  
    2323                                        error('constructor not supported');
    2424                        end
     25                end % }}}
     26                function disp(self) % {{{
     27                        disp(sprintf('   timestepping parameters:'));
     28
     29                        unit = 'yr';
     30                        fielddisplay(self,'start_time',['simulation starting time [' unit ']']);
     31                        fielddisplay(self,'final_time',['final time to stop the simulation [' unit ']']);
     32                        fielddisplay(self,'time_step',['length of time steps [' unit ']']);
     33                        fielddisplay(self,'interp_forcing','interpolate in time between requested forcing values? (0 or 1)');
     34                        fielddisplay(self,'cycle_forcing','cycle through forcing? (0 or 1)');
     35                        fielddisplay(self,'coupling_time',['length of coupling time step with ocean model [' unit ']']);
     36
    2537                end % }}}
    2638                function self = setdefaultparameters(self) % {{{
     
    5163                        end
    5264                end % }}}
    53                 function disp(self) % {{{
    54                         disp(sprintf('   timestepping parameters:'));
    55 
    56                         unit = 'yr';
    57                         fielddisplay(self,'start_time',['simulation starting time [' unit ']']);
    58                         fielddisplay(self,'final_time',['final time to stop the simulation [' unit ']']);
    59                         fielddisplay(self,'time_step',['length of time steps [' unit ']']);
    60                         fielddisplay(self,'interp_forcing','interpolate in time between requested forcing values ? (0 or 1)');
    61                         fielddisplay(self,'cycle_forcing','cycle through forcing ? (0 or 1)');
    62                         fielddisplay(self,'coupling_time',['length of coupling time step with ocean model  [' unit ']']);
    63 
    64                 end % }}}
    6565                function marshall(self,prefix,md,fid) % {{{
    6666
  • issm/trunk-jpl/src/m/classes/timestepping.py

    r26208 r26303  
    55
    66class timestepping(object):
    7     """
    8     TIMESTEPPING Class definition
     7    """TIMESTEPPING Class definition
    98
    10        Usage:
    11           timestepping = timestepping()
     9    Usage:
     10        timestepping = timestepping()
    1211    """
    1312
    14     def __init__(self):  # {{{
    15         self.start_time = 0.
    16         self.final_time = 0.
    17         self.time_step = 0.
     13    def __init__(self, *args): #{{{
     14        self.start_time = 0
     15        self.final_time = 0
     16        self.time_step = 0
    1817        self.interp_forcing = 1
    1918        self.cycle_forcing = 0
    20         self.coupling_time = 0.
     19        self.coupling_time = 0
    2120
    22     #set defaults
    23         self.setdefaultparameters()
    24 
     21        if len(args) == 0:
     22            self.setdefaultparameters()
     23        else:
     24            raise RuntimeError('constructor not supported')
    2525    #}}}
    2626
    27     def __repr__(self):  # {{{
    28         string = "   timestepping parameters:"
    29         string = "%s\n%s" % (string, fielddisplay(self, "start_time", "simulation starting time [yr]"))
    30         string = "%s\n%s" % (string, fielddisplay(self, "final_time", "final time to stop the simulation [yr]"))
    31         string = "%s\n%s" % (string, fielddisplay(self, "time_step", "length of time steps [yr]"))
    32         string = "%s\n%s" % (string, fielddisplay(self, "interp_forcing", "interpolate in time between requested forcing values ? (0 or 1)"))
    33         string = "%s\n%s" % (string, fielddisplay(self, "cycle_forcing", "cycle through forcing ? (0 or 1)"))
    34         string = "%s\n%s" % (string, fielddisplay(self, "coupling_time", "length of coupling time steps with ocean model [yr]"))
    35         return string
     27    def __repr__(self): #{{{
     28        s = '   timestepping parameters:\n'
     29        unit = 'yr'
     30        s += '{}\n'.format(fielddisplay(self, 'start_time', 'simulation starting time [' + unit + ']'))
     31        s += '{}\n'.format(fielddisplay(self, 'final_time', 'final time to stop the simulation [' + unit + ']'))
     32        s += '{}\n'.format(fielddisplay(self, 'time_step', 'length of time steps [' + unit + ']'))
     33        s += '{}\n'.format(fielddisplay(self, 'interp_forcing', 'interpolate in time between requested forcing values? (0 or 1)'))
     34        s += '{}\n'.format(fielddisplay(self, 'cycle_forcing', 'cycle through forcing? (0 or 1)'))
     35        s += '{}\n'.format(fielddisplay(self, 'coupling_time', 'length of coupling time steps with ocean model [' + unit + ']'))
     36        return s
    3637    #}}}
    3738
    38     def setdefaultparameters(self):  # {{{
    39         #time between 2 time steps
    40         self.time_step = 1. / 2.
    41         #final time
    42         self.final_time = 10. * self.time_step
    43         #should we interpolate forcing between timesteps?
     39    def setdefaultparameters(self): #{{{
     40        # Time between 2 time steps
     41        self.time_step = 1 / 2
     42
     43        # Final time
     44        self.final_time = 10 * self.time_step
     45
     46        # Should we interpolate forcing between timesteps?
    4447        self.interp_forcing = 1
    4548        self.cycle_forcing = 0
     
    4851    #}}}
    4952
    50     def checkconsistency(self, md, solution, analyses):  # {{{
    51 
     53    def checkconsistency(self, md, solution, analyses): #{{{
    5254        md = checkfield(md, 'fieldname', 'timestepping.start_time', 'numel', [1], 'NaN', 1, 'Inf', 1)
    5355        md = checkfield(md, 'fieldname', 'timestepping.final_time', 'numel', [1], 'NaN', 1, 'Inf', 1)
    5456        md = checkfield(md, 'fieldname', 'timestepping.time_step', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1)
    55         if self.final_time - self.start_time < 0:
    56             md.checkmessage("timestepping.final_time should be larger than timestepping.start_time")
    57             md = checkfield(md, 'fieldname', 'timestepping.coupling_time', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1)
    58         md = checkfield(md, 'fieldname', 'timestepping.interp_forcing', 'numel', [1], 'values', [0, 1])
    5957        md = checkfield(md, 'fieldname', 'timestepping.cycle_forcing', 'numel', [1], 'values', [0, 1])
     58        if (self.final_time - self.start_time) < 0:
     59            md.checkmessage('timestepping.final_time should be larger than timestepping.start_time')
     60        if solution == 'TransientSolution':
     61            md = checkfield(md, 'fieldname', 'timestepping.time_step', 'numel', [1], '>', 0, 'NaN', 1, 'Inf', 1)
     62            md = checkfield(md, 'fieldname', 'timestepping.time_step', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1)
    6063
    6164        return md
    62     # }}}
     65    #}}}
    6366
    64     def marshall(self, prefix, md, fid):  # {{{
    65 
    66         yts = md.constants.yts
     67    def marshall(self, prefix, md, fid): #{{{
     68        scale = md.constants.yts
    6769        WriteData(fid, prefix, 'name', 'md.timestepping.type', 'data', 1, 'format', 'Integer')
    68         WriteData(fid, prefix, 'object', self, 'fieldname', 'start_time', 'format', 'Double', 'scale', yts)
    69         WriteData(fid, prefix, 'object', self, 'fieldname', 'final_time', 'format', 'Double', 'scale', yts)
    70         WriteData(fid, prefix, 'object', self, 'fieldname', 'time_step', 'format', 'Double', 'scale', yts)
     70        WriteData(fid, prefix, 'object', self, 'fieldname', 'start_time', 'format', 'Double', 'scale', scale)
     71        WriteData(fid, prefix, 'object', self, 'fieldname', 'final_time', 'format', 'Double', 'scale', scale)
     72        WriteData(fid, prefix, 'object', self, 'fieldname', 'time_step', 'format', 'Double', 'scale', scale)
    7173        WriteData(fid, prefix, 'object', self, 'fieldname', 'interp_forcing', 'format', 'Boolean')
    7274        WriteData(fid, prefix, 'object', self, 'fieldname', 'cycle_forcing', 'format', 'Boolean')
    73         WriteData(fid, prefix, 'object', self, 'fieldname', 'coupling_time', 'format', 'Double', 'scale', yts)
    74     # }}}
     75        WriteData(fid, prefix, 'object', self, 'fieldname', 'coupling_time', 'format', 'Double', 'scale', scale)
     76    #}}}
  • issm/trunk-jpl/test/NightlyRun/test2002.m

    r26301 r26303  
    9898Seustatic=md.results.TransientSolution.Sealevel;
    9999Beustatic=md.results.TransientSolution.Bed;
    100 pause
    101100
    102101%eustatic + selfattraction run:
  • issm/trunk-jpl/test/NightlyRun/test2002.py

    r26301 r26303  
    2323
    2424# Geometry for the bed, arbitrary thickness of 100
    25 md.geometry.bed = np.zeros((md.mesh.numberofvertices, 1))
     25md.geometry.bed = np.zeros((md.mesh.numberofvertices, ))
    2626md.geometry.base = md.geometry.bed
    27 md.geometry.thickness = 100 * np.ones((md.mesh.numberofvertices, 1))
     27md.geometry.thickness = 100 * np.ones((md.mesh.numberofvertices, ))
    2828md.geometry.surface = md.geometry.bed + md.geometry.thickness
    2929
     
    5555# Mask: {{{
    5656mask = gmtmask(md.mesh.lat, md.mesh.long)
    57 oceanmask = -1 * np.ones((md.mesh.numberofvertices, 1))
     57oceanmask = -1 * np.ones((md.mesh.numberofvertices, ))
    5858pos = np.where(mask == 0)[0]
    5959oceanmask[pos] = 1
    6060
    61 icemask = np.ones((md.mesh.numberofvertices, 1))
     61icemask = np.ones((md.mesh.numberofvertices, ))
    6262icemask[md.mesh.elements[posant]] = -1
    6363icemask[md.mesh.elements[posgre]] = -1
     
    7474
    7575# Masstransport
    76 md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices, 1))
    77 md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices, 1))
    78 md.initialization.vx = np.zeros((md.mesh.numberofvertices, 1))
    79 md.initialization.vy = np.zeros((md.mesh.numberofvertices, 1))
    80 md.initialization.sealevel = np.zeros((md.mesh.numberofvertices, 1))
     76md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices, ))
     77md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices, ))
     78md.initialization.vx = np.zeros((md.mesh.numberofvertices, ))
     79md.initialization.vy = np.zeros((md.mesh.numberofvertices, ))
     80md.initialization.sealevel = np.zeros((md.mesh.numberofvertices, ))
    8181md.initialization.str = 0
    8282
  • issm/trunk-jpl/test/NightlyRun/test2003.m

    r26296 r26303  
    1212md.geometry.surface=md.geometry.bed+md.geometry.thickness;
    1313
    14 
    1514%parameterize slc solution:
    1615%solidearth loading:  {{{
    1716md.masstransport.spcthickness=[md.geometry.thickness;0];
    1817md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
    19 
    2018
    2119xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3;
     
    4442md.mask.ocean_levelset=oceanmask;
    4543
    46 % use model representation of ocen area (not the true area)
     44% use model representation of ocean area (not the true area)
    4745md.solidearth.settings.ocean_area_scaling = 0;
    4846
  • issm/trunk-jpl/test/NightlyRun/test2003.py

    r25956 r26303  
    1212
    1313
    14 #mesh earth:
     14# Mesh earth
     15#
     16# NOTE: In MATLAB, we currently use cached mesh to account for differences in
     17# mesh generated under Linux versus under macOS
     18#
    1519md = model()
    16 md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 1000.)  #1000 km resolution mesh
     20md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 700.) # 700 km resolution mesh
    1721
    18 #parameterize solidearth solution:
     22# Geometry for the bed, arbitrary thickness of 100
     23md.geometry.bed = -1 * np.ones(md.mesh.numberofvertices)
     24md.geometry.base = md.geometry.bed
     25md.geometry.thickness = 1000 * np.ones(md.mesh.numberofvertices)
     26md.geometry.surface = md.geometry.bed + md.geometry.thickness
     27
     28# Parameterize slc solution:
    1929#solidearth loading:  {{{
    20 md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, ))
    21 md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, ))
    22 md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, ))
    23 md.dsl.sea_surface_height_change_above_geoid=np.zeros((md.mesh.numberofvertices+1, ))
    24 md.dsl.sea_water_pressure_change_at_sea_floor=np.zeros((md.mesh.numberofvertices+1, ))
     30md.masstransport.spcthickness = np.append(md.geometry.thickness, 0)
     31md.smb.mass_balance = np.zeros(md.mesh.numberofvertices)
    2532
    26 #antarctica
    27 #Access every element in lat using the indices in elements
    28 # - 1 to convert to base 0 indexing, 1 (not 2, in matlab) to sum over rows
    29 late = md.mesh.lat[md.mesh.elements - 1].sum(axis=1)/ 3
    30 longe = md.mesh.long[md.mesh.elements - 1].sum(axis=1)/ 3
    31 pos = np.intersect1d(np.array(np.where(late < -75)[0]), np.array(np.where(longe < 0)[0]))
    32 md.solidearth.surfaceload.icethicknesschange[pos] = -1
     33xe = md.mesh.x[md.mesh.elements - 1].sum(axis=1) / 3
     34ye = md.mesh.y[md.mesh.elements - 1].sum(axis=1) / 3
     35ze = md.mesh.z[md.mesh.elements - 1].sum(axis=1) / 3
     36re = pow((pow(xe, 2) + pow(ye, 2) + pow(ze, 2)), 0.5)
    3337
    34 #elastic loading from love numbers:
    35 md.solidearth.lovenumbers = lovenumbers('maxdeg', 1000)
     38late = asind(ze / re)
     39longe = atan2d(ye, xe)
     40# Greenland
     41pos = np.where(np.logical_and.reduce((late > 60, late < 90, longe > -75, longe < -15)))[0]
     42md.masstransport.spcthickness[md.mesh.elements[pos]] = md.masstransport.spcthickness[md.mesh.elements[pos]] - 1000
     43posice = pos
     44
     45# Elastic loading from love numbers
     46md.solidearth.lovenumbers = lovenumbers('maxdeg', 100)
    3647#}}}
    3748
    38 #mask: {{{
     49# Mask: {{{
    3950mask = gmtmask(md.mesh.lat, md.mesh.long)
    4051icemask = np.ones(md.mesh.numberofvertices)
     52icemask[md.mesh.elements[posice]] = -1
     53md.mask.ice_levelset = icemask
     54oceanmask = -1 * np.ones(md.mesh.numberofvertices)
    4155pos = np.where(mask == 0)[0]
    42 icemask[pos] = -1
    43 pos = np.where(mask[md.mesh.elements - 1].sum(axis=1) < 3)[0]
    44 icemask[md.mesh.elements[pos, :] - 1] = -1
    45 md.mask.ice_levelset = icemask
    46 md.mask.ocean_levelset = -icemask
     56oceanmask[pos] = 1
     57md.mask.ocean_levelset = oceanmask
    4758
    48 #make sure that the elements that have loads are fully grounded
    49 pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0]
    50 md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1
    51 
    52 #make sure wherever there is an ice load, that the mask is set to ice:
    53 #pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] # Do we need to do this twice?
    54 md.mask.ice_levelset[md.mesh.elements[pos, :] - 1] = -1
    55 # }}}
    56 
    57 # use model representation of ocen area (not the true area)
     59# Use model representation of ocean area (not the true area)
    5860md.solidearth.settings.ocean_area_scaling = 0
    5961
    60 #geometry
    61 di = md.materials.rho_ice / md.materials.rho_water
    62 md.geometry.thickness = np.ones(md.mesh.numberofvertices)
    63 md.geometry.surface = (1 - di) * np.zeros(md.mesh.numberofvertices)
    64 md.geometry.base = md.geometry.surface - md.geometry.thickness
    65 md.geometry.bed = md.geometry.base
     62# Materials
     63md.initialization.temperature = 273.25 * np.ones(md.mesh.numberofvertices)
     64md.initialization.sealevel = np.zeros(md.mesh.numberofvertices)
     65md.initialization.str = 0
    6666
    67 #materials
    68 md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices, ))
    69 md.materials.rheology_B = paterson(md.initialization.temperature)
    70 md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, ))
     67md.basalforcings.groundedice_melting_rate = np.zeros(md.mesh.numberofvertices)
     68md.basalforcings.floatingice_melting_rate = np.zeros(md.mesh.numberofvertices)
     69md.initialization.vx = np.zeros(md.mesh.numberofvertices)
     70md.initialization.vy = np.zeros(md.mesh.numberofvertices)
    7171
    72 #Miscellaneous
     72# Miscellaneous
    7373md.miscellaneous.name = 'test2003'
    7474
    75 #Solution parameters
     75# Solution parameters
    7676md.solidearth.settings.reltol = np.nan
    7777md.solidearth.settings.abstol = 1e-3
    78 md.solidearth.settings.computesealevelchange = 1
     78md.solidearth.settings.sealevelloading = 0
     79md.solidearth.settings.grdocean = 0
     80md.solidearth.settings.isgrd = 1
     81md.solidearth.settings.ocean_area_scaling = 0
     82md.solidearth.settings.grdmodel = 1
     83md.solidearth.settings.horiz = 1
     84md.solidearth.requested_outputs = ['Sealevel', 'Bed', 'BedEast', 'BedNorth']
    7985
    80 #eustatic + rigid + elastic run:
    81 md.solidearth.settings.rigid = 1
     86# Physics
     87md.transient.issmb = 0
     88md.transient.isstressbalance = 0
     89md.transient.isthermal = 0
     90md.transient.ismasstransport = 1
     91md.transient.isslc = 1
     92
     93md.timestepping.start_time = 0
     94md.timestepping.time_step = 1
     95md.timestepping.final_time = 1
     96
     97# Eustatic + selfattraction + elastic run:
     98md.solidearth.settings.selfattraction = 1
    8299md.solidearth.settings.elastic = 1
    83100md.solidearth.settings.rotation = 0
     101md.solidearth.settings.viscous = 0
    84102md.cluster = generic('name', gethostname(), 'np', 3)
    85103#md.verbose = verbose('111111111')
    86 md = solve(md, 'Sealevelrise')
    87 SnoRotation = md.results.SealevelriseSolution.Sealevel
     104md = solve(md, 'Transient')
     105SnoRotation = md.results.TransientSolution.Sealevel
     106BUnoRotation = md.results.TransientSolution.Bed
     107BEnoRotation = md.results.TransientSolution.BedEast
     108BNnoRotation = md.results.TransientSolution.BedNorth
    88109
    89 #eustatic + rigid + elastic + rotation run:
    90 md.solidearth.settings.rigid = 1
     110# Eustatic + selfattraction + elastic + rotation run
     111md.solidearth.settings.selfattraction = 1
    91112md.solidearth.settings.elastic = 1
    92113md.solidearth.settings.rotation = 1
     114md.solidearth.settings.viscous = 0
    93115md.cluster = generic('name', gethostname(), 'np', 3)
    94116#md.verbose = verbose('111111111')
    95 md = solve(md, 'Sealevelrise')
    96 SRotation = md.results.SealevelriseSolution.Sealevel
     117md = solve(md, 'Transient')
     118SRotation = md.results.TransientSolution.Sealevel
     119BURotation = md.results.TransientSolution.Bed
     120BERotation = md.results.TransientSolution.BedEast
     121BNRotation = md.results.TransientSolution.BedNorth
    97122
    98 #Fields and tolerances to track changes
    99 field_names = ['noRotation', 'Rotation']
    100 field_tolerances = [1e-13, 1e-13]
    101 field_values = [SnoRotation, SRotation]
     123# Fields and tolerances to track changes
     124field_names = ['Sealevel', 'Uplift', 'NorthDisplacement', 'EastDisplacement']
     125field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13]
     126field_values = [SRotation - SnoRotation, BURotation - BUnoRotation, BNRotation - BNnoRotation,BERotation - BEnoRotation]
Note: See TracChangeset for help on using the changeset viewer.