Changeset 26303
- Timestamp:
- 06/08/21 00:40:42 (4 years ago)
- Location:
- issm/trunk-jpl
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/classes/timestepping.m
r26208 r26303 23 23 error('constructor not supported'); 24 24 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 25 37 end % }}} 26 38 function self = setdefaultparameters(self) % {{{ … … 51 63 end 52 64 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 % }}}65 65 function marshall(self,prefix,md,fid) % {{{ 66 66 -
issm/trunk-jpl/src/m/classes/timestepping.py
r26208 r26303 5 5 6 6 class timestepping(object): 7 """ 8 TIMESTEPPING Class definition 7 """TIMESTEPPING Class definition 9 8 10 11 9 Usage: 10 timestepping = timestepping() 12 11 """ 13 12 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 18 17 self.interp_forcing = 1 19 18 self.cycle_forcing = 0 20 self.coupling_time = 0 .19 self.coupling_time = 0 21 20 22 #set defaults 23 self.setdefaultparameters() 24 21 if len(args) == 0: 22 self.setdefaultparameters() 23 else: 24 raise RuntimeError('constructor not supported') 25 25 #}}} 26 26 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 36 37 #}}} 37 38 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? 44 47 self.interp_forcing = 1 45 48 self.cycle_forcing = 0 … … 48 51 #}}} 49 52 50 def checkconsistency(self, md, solution, analyses): # {{{ 51 53 def checkconsistency(self, md, solution, analyses): #{{{ 52 54 md = checkfield(md, 'fieldname', 'timestepping.start_time', 'numel', [1], 'NaN', 1, 'Inf', 1) 53 55 md = checkfield(md, 'fieldname', 'timestepping.final_time', 'numel', [1], 'NaN', 1, 'Inf', 1) 54 56 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])59 57 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) 60 63 61 64 return md 62 # 65 #}}} 63 66 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 67 69 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) 71 73 WriteData(fid, prefix, 'object', self, 'fieldname', 'interp_forcing', 'format', 'Boolean') 72 74 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 98 98 Seustatic=md.results.TransientSolution.Sealevel; 99 99 Beustatic=md.results.TransientSolution.Bed; 100 pause101 100 102 101 %eustatic + selfattraction run: -
issm/trunk-jpl/test/NightlyRun/test2002.py
r26301 r26303 23 23 24 24 # Geometry for the bed, arbitrary thickness of 100 25 md.geometry.bed = np.zeros((md.mesh.numberofvertices, 1))25 md.geometry.bed = np.zeros((md.mesh.numberofvertices, )) 26 26 md.geometry.base = md.geometry.bed 27 md.geometry.thickness = 100 * np.ones((md.mesh.numberofvertices, 1))27 md.geometry.thickness = 100 * np.ones((md.mesh.numberofvertices, )) 28 28 md.geometry.surface = md.geometry.bed + md.geometry.thickness 29 29 … … 55 55 # Mask: {{{ 56 56 mask = gmtmask(md.mesh.lat, md.mesh.long) 57 oceanmask = -1 * np.ones((md.mesh.numberofvertices, 1))57 oceanmask = -1 * np.ones((md.mesh.numberofvertices, )) 58 58 pos = np.where(mask == 0)[0] 59 59 oceanmask[pos] = 1 60 60 61 icemask = np.ones((md.mesh.numberofvertices, 1))61 icemask = np.ones((md.mesh.numberofvertices, )) 62 62 icemask[md.mesh.elements[posant]] = -1 63 63 icemask[md.mesh.elements[posgre]] = -1 … … 74 74 75 75 # 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))76 md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices, )) 77 md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices, )) 78 md.initialization.vx = np.zeros((md.mesh.numberofvertices, )) 79 md.initialization.vy = np.zeros((md.mesh.numberofvertices, )) 80 md.initialization.sealevel = np.zeros((md.mesh.numberofvertices, )) 81 81 md.initialization.str = 0 82 82 -
issm/trunk-jpl/test/NightlyRun/test2003.m
r26296 r26303 12 12 md.geometry.surface=md.geometry.bed+md.geometry.thickness; 13 13 14 15 14 %parameterize slc solution: 16 15 %solidearth loading: {{{ 17 16 md.masstransport.spcthickness=[md.geometry.thickness;0]; 18 17 md.smb.mass_balance=zeros(md.mesh.numberofvertices,1); 19 20 18 21 19 xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3; … … 44 42 md.mask.ocean_levelset=oceanmask; 45 43 46 % use model representation of oce n area (not the true area)44 % use model representation of ocean area (not the true area) 47 45 md.solidearth.settings.ocean_area_scaling = 0; 48 46 -
issm/trunk-jpl/test/NightlyRun/test2003.py
r25956 r26303 12 12 13 13 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 # 15 19 md = model() 16 md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 1000.) #1000 km resolution mesh20 md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 700.) # 700 km resolution mesh 17 21 18 #parameterize solidearth solution: 22 # Geometry for the bed, arbitrary thickness of 100 23 md.geometry.bed = -1 * np.ones(md.mesh.numberofvertices) 24 md.geometry.base = md.geometry.bed 25 md.geometry.thickness = 1000 * np.ones(md.mesh.numberofvertices) 26 md.geometry.surface = md.geometry.bed + md.geometry.thickness 27 28 # Parameterize slc solution: 19 29 #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, )) 30 md.masstransport.spcthickness = np.append(md.geometry.thickness, 0) 31 md.smb.mass_balance = np.zeros(md.mesh.numberofvertices) 25 32 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 33 xe = md.mesh.x[md.mesh.elements - 1].sum(axis=1) / 3 34 ye = md.mesh.y[md.mesh.elements - 1].sum(axis=1) / 3 35 ze = md.mesh.z[md.mesh.elements - 1].sum(axis=1) / 3 36 re = pow((pow(xe, 2) + pow(ye, 2) + pow(ze, 2)), 0.5) 33 37 34 #elastic loading from love numbers: 35 md.solidearth.lovenumbers = lovenumbers('maxdeg', 1000) 38 late = asind(ze / re) 39 longe = atan2d(ye, xe) 40 # Greenland 41 pos = np.where(np.logical_and.reduce((late > 60, late < 90, longe > -75, longe < -15)))[0] 42 md.masstransport.spcthickness[md.mesh.elements[pos]] = md.masstransport.spcthickness[md.mesh.elements[pos]] - 1000 43 posice = pos 44 45 # Elastic loading from love numbers 46 md.solidearth.lovenumbers = lovenumbers('maxdeg', 100) 36 47 #}}} 37 48 38 # mask:{{{49 # Mask: {{{ 39 50 mask = gmtmask(md.mesh.lat, md.mesh.long) 40 51 icemask = np.ones(md.mesh.numberofvertices) 52 icemask[md.mesh.elements[posice]] = -1 53 md.mask.ice_levelset = icemask 54 oceanmask = -1 * np.ones(md.mesh.numberofvertices) 41 55 pos = 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 56 oceanmask[pos] = 1 57 md.mask.ocean_levelset = oceanmask 47 58 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) 58 60 md.solidearth.settings.ocean_area_scaling = 0 59 61 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 63 md.initialization.temperature = 273.25 * np.ones(md.mesh.numberofvertices) 64 md.initialization.sealevel = np.zeros(md.mesh.numberofvertices) 65 md.initialization.str = 0 66 66 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, ))67 md.basalforcings.groundedice_melting_rate = np.zeros(md.mesh.numberofvertices) 68 md.basalforcings.floatingice_melting_rate = np.zeros(md.mesh.numberofvertices) 69 md.initialization.vx = np.zeros(md.mesh.numberofvertices) 70 md.initialization.vy = np.zeros(md.mesh.numberofvertices) 71 71 72 # Miscellaneous72 # Miscellaneous 73 73 md.miscellaneous.name = 'test2003' 74 74 75 # Solution parameters75 # Solution parameters 76 76 md.solidearth.settings.reltol = np.nan 77 77 md.solidearth.settings.abstol = 1e-3 78 md.solidearth.settings.computesealevelchange = 1 78 md.solidearth.settings.sealevelloading = 0 79 md.solidearth.settings.grdocean = 0 80 md.solidearth.settings.isgrd = 1 81 md.solidearth.settings.ocean_area_scaling = 0 82 md.solidearth.settings.grdmodel = 1 83 md.solidearth.settings.horiz = 1 84 md.solidearth.requested_outputs = ['Sealevel', 'Bed', 'BedEast', 'BedNorth'] 79 85 80 #eustatic + rigid + elastic run: 81 md.solidearth.settings.rigid = 1 86 # Physics 87 md.transient.issmb = 0 88 md.transient.isstressbalance = 0 89 md.transient.isthermal = 0 90 md.transient.ismasstransport = 1 91 md.transient.isslc = 1 92 93 md.timestepping.start_time = 0 94 md.timestepping.time_step = 1 95 md.timestepping.final_time = 1 96 97 # Eustatic + selfattraction + elastic run: 98 md.solidearth.settings.selfattraction = 1 82 99 md.solidearth.settings.elastic = 1 83 100 md.solidearth.settings.rotation = 0 101 md.solidearth.settings.viscous = 0 84 102 md.cluster = generic('name', gethostname(), 'np', 3) 85 103 #md.verbose = verbose('111111111') 86 md = solve(md, 'Sealevelrise') 87 SnoRotation = md.results.SealevelriseSolution.Sealevel 104 md = solve(md, 'Transient') 105 SnoRotation = md.results.TransientSolution.Sealevel 106 BUnoRotation = md.results.TransientSolution.Bed 107 BEnoRotation = md.results.TransientSolution.BedEast 108 BNnoRotation = md.results.TransientSolution.BedNorth 88 109 89 # eustatic + rigid + elastic + rotation run:90 md.solidearth.settings. rigid= 1110 # Eustatic + selfattraction + elastic + rotation run 111 md.solidearth.settings.selfattraction = 1 91 112 md.solidearth.settings.elastic = 1 92 113 md.solidearth.settings.rotation = 1 114 md.solidearth.settings.viscous = 0 93 115 md.cluster = generic('name', gethostname(), 'np', 3) 94 116 #md.verbose = verbose('111111111') 95 md = solve(md, 'Sealevelrise') 96 SRotation = md.results.SealevelriseSolution.Sealevel 117 md = solve(md, 'Transient') 118 SRotation = md.results.TransientSolution.Sealevel 119 BURotation = md.results.TransientSolution.Bed 120 BERotation = md.results.TransientSolution.BedEast 121 BNRotation = md.results.TransientSolution.BedNorth 97 122 98 # Fields and tolerances to track changes99 field_names = [' noRotation', 'Rotation']100 field_tolerances = [1e-13, 1e-13 ]101 field_values = [S noRotation, SRotation]123 # Fields and tolerances to track changes 124 field_names = ['Sealevel', 'Uplift', 'NorthDisplacement', 'EastDisplacement'] 125 field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13] 126 field_values = [SRotation - SnoRotation, BURotation - BUnoRotation, BNRotation - BNnoRotation,BERotation - BEnoRotation]
Note:
See TracChangeset
for help on using the changeset viewer.