PKrCHXI{wwstagpy/commands.py"""definition of each subcommands""" from . import constants, misc, field, rprof, time_series, plates def field_cmd(args): """plot snapshots of fields""" misc.parse_timesteps(args) misc.plot_backend(args) if args.plot is not None: for var, meta in constants.FIELD_VAR_LIST.items(): misc.set_arg(args, meta.arg, var in args.plot) field.field_cmd(args) def rprof_cmd(args): """plot radial profiles""" misc.parse_timesteps(args) misc.plot_backend(args) if args.plot is not None: for var, meta in constants.RPROF_VAR_LIST.items(): if var in args.plot: misc.set_arg(args, meta.arg, True) misc.set_arg(args, meta.min_max, False) else: misc.set_arg(args, meta.arg, False) if meta.min_max and var.upper() in args.plot: misc.set_arg(args, meta.arg, True) misc.set_arg(args, meta.min_max, True) rprof.rprof_cmd(args) def time_cmd(args): """plot time series""" misc.parse_timesteps(args) if args.compstat: args.matplotback = None misc.plot_backend(args) time_series.time_cmd(args) def plates_cmd(args): """plate analysis""" misc.parse_timesteps(args) misc.plot_backend(args) plates.plates_cmd(args) def var_cmd(_): """display a list of available variables""" print('field:') print(*('{}: {}'.format(v, m.name) for v, m in constants.FIELD_VAR_LIST.items()), sep='\n') print() print('rprof:') print(*('{}: {}'.format(v, m.name) for v, m in constants.RPROF_VAR_LIST.items()), sep='\n') PKmFHHp= stagpy/time_series.py"""Plots time series of temperature and heat fluxes outputs from stagyy. Author: Stephane Labrosse with inputs from Martina Ulvrova and Adrien Morison Date: 2015/11/27 """ import numpy as np from math import sqrt from .stagdata import TimeData def find_nearest(array, value): """Find the data point nearest to value""" idx = (np.abs(array - value)).argmin() return array[idx] def time_cmd(args): """plot temporal series""" plt = args.plt lwdth = args.linewidth ftsz = args.fontsize rab = args.par_nml['refstate']['Ra0'] rah = args.par_nml['refstate']['Rh'] botpphase = args.par_nml['boundaries']['BotPphase'] time_data = TimeData(args) colnames, data = time_data.colnames, time_data.data ntot = len(data) spherical = args.par_nml['geometry']['shape'].lower() == 'spherical' if spherical: rcmb = args.par_nml['geometry']['r_cmb'] rmin = rcmb rmax = rmin + 1 coefb = 1 # rb**2*4*pi coefs = (rmax / rmin)**2 # *4*pi volume = rmin * (1 - (rmax / rmin)**3) / 3 # *4*pi/3 else: rcmb = 0. coefb = 1. coefs = 1. volume = 1. time = data[:, 1] dtdt = (data[2:ntot - 1, 5] - data[0:ntot - 3, 5]) /\ (data[2:ntot - 1, 1] - data[0:ntot - 3, 1]) ebalance = data[1:ntot - 2, 2] * coefs - data[1:ntot - 2, 3] * coefb\ - volume * dtdt fig = plt.figure(figsize=(30, 10)) plt.subplot(2, 1, 1) plt.plot(time, data[:, 2] * coefs, 'b', label='Surface', linewidth=lwdth) plt.plot(time, data[:, 3] * coefb, 'r', label='Bottom', linewidth=lwdth) plt.plot(time[1:ntot - 2:], ebalance, 'g', label='Energy balance', linewidth=lwdth) plt.ylabel('Heat flow', fontsize=ftsz) plt.legend = plt.legend(loc='upper right', shadow=False, fontsize=ftsz) plt.legend.get_frame().set_facecolor('white') plt.xticks(fontsize=ftsz) plt.yticks(fontsize=ftsz) if args.annottmin: plt.annotate('tminT', xy=(args.tmint, 0), xytext=(args.tmint, -10), arrowprops={'facecolor': 'black'}) plt.annotate('tminC', xy=(args.tminc, 0), xytext=(args.tminc, 10), arrowprops={'facecolor': 'black'}) plt.subplot(2, 1, 2) plt.plot(time, data[:, 5], 'k', linewidth=lwdth) plt.xlabel('Time', fontsize=ftsz) plt.ylabel('Mean temperature', fontsize=ftsz) plt.xticks(fontsize=ftsz) plt.yticks(fontsize=ftsz) plt.savefig("flux_time.pdf", format='PDF') if not args.compstat: return None coords = [] print('right click to select starting time of statistics computations') def onclick(event): """get position and button from mouse click""" ixc, iyc = event.xdata, event.ydata button = event.button # assign global variable to access outside of function if button == 3: coords.append((ixc, iyc)) # Disconnect after 1 clicks if len(coords) == 1: fig.canvas.mpl_disconnect(cid) plt.close(1) return # Call click func cid = fig.canvas.mpl_connect('button_press_event', onclick) plt.show() moy = [] rms = [] ebal = [] rms_ebal = [] ch1 = np.where(time == (find_nearest(time, coords[0][0]))) print('Statistics computed from t =' + str(time[ch1[0][0]])) for num in range(2, len(colnames)): moy.append(np.trapz(data[ch1[0][0]:ntot - 1, num], x=time[ch1[0][0]:ntot - 1]) / (time[ntot - 1] - time[ch1[0][0]])) rms.append(sqrt(np.trapz((data[ch1[0][0]:ntot - 1, num] - moy[num - 2])**2, x=time[ch1[0][0]:ntot - 1]) / (time[ntot - 1] - time[ch1[0][0]]))) print(colnames[num], '=', moy[num - 2], 'pm', rms[num - 2]) ebal.append(np.trapz(ebalance[ch1[0][0] - 1:ntot - 3], x=time[ch1[0][0]:ntot - 2]) / (time[ntot - 2] - time[ch1[0][0]])) rms_ebal.append(sqrt(np.trapz( (ebalance[ch1[0][0] - 1:ntot - 3] - ebal)**2, x=time[ch1[0][0]:ntot - 2]) / (time[ntot - 2] - time[ch1[0][0]]))) print('Energy balance', ebal, 'pm', rms_ebal) results = moy + ebal + rms + rms_ebal with open('Stats.dat', 'w') as out_file: out_file.write("%10.5e %10.5e %10.5e " % (rab, rah, botpphase)) for item in results: out_file.write("%10.5e " % item) out_file.write("\n") PK4HWWstagpy/__main__.py# PYTHON_ARGCOMPLETE_OK """Makes stagpy/ callable""" from .stagpy import main main() PKXILHz3z3stagpy/config.py"""Handle configuration of StagPy Create the cmd line argument parser and deal with the config file """ from collections import OrderedDict, namedtuple from os import mkdir from subprocess import call import argcomplete import argparse import configparser import os.path import shlex from . import commands, parfile from .constants import CONFIG_DIR CONFIG_FILE = os.path.join(CONFIG_DIR, 'config') Conf = namedtuple('ConfigEntry', ['default', 'cmd_arg', 'shortname', 'kwargs', 'conf_arg', 'help_string']) CORE = OrderedDict(( ('path', Conf('./', True, 'p', {}, True, 'StagYY run directory')), ('name', Conf('test', True, None, {}, True, 'StagYY generic output file name')), ('outname', Conf('stagpy', True, 'n', {}, True, 'StagPy generic output file name')), ('geometry', Conf('annulus', True, 'g', {'choices': ['annulus']}, True, 'geometry of the domain')), ('timestep', Conf('100', True, 's', {}, True, 'timestep slice')), ('xkcd', Conf(False, True, None, {}, True, 'use the xkcd style')), ('pdf', Conf(False, True, None, {}, True, 'produce non-rasterized pdf (slow!)')), ('fontsize', Conf(16, False, None, {}, True, 'font size')), ('linewidth', Conf(2, False, None, {}, True, 'line width')), ('matplotback', Conf('agg', False, None, {}, True, 'graphical backend')), ('useseaborn', Conf(True, False, None, {}, True, 'use or not seaborn')), )) FIELD = OrderedDict(( ('plot', Conf(None, True, 'o', {'nargs': '?', 'const': '', 'type': str}, False, 'specify which variable to plot')), ('plot_temperature', Conf(True, False, None, {}, True, 'temperature scalar field')), ('plot_xvelo', Conf(False, False, None, {}, True, 'x velocity scalar field')), ('plot_yvelo', Conf(False, False, None, {}, True, 'y velocity scalar field')), ('plot_zvelo', Conf(False, False, None, {}, True, 'z velocity scalar field')), ('plot_pressure', Conf(True, False, None, {}, True, 'pressure scalar field')), ('plot_stream', Conf(True, False, None, {}, True, 'stream function scalar field')), ('plot_composition', Conf(False, False, None, {}, True, 'composition scalar field')), ('plot_viscosity', Conf(False, False, None, {}, True, 'viscosity scalar field')), ('plot_density', Conf(False, False, None, {}, True, 'density scalar field')), ('plot_water', Conf(False, False, None, {}, True, 'water concentration scalar field')), ('plot_age', Conf(False, False, None, {}, True, 'age scalar field')), ('shrinkcb', Conf(0.5, False, None, {}, True, 'color bar shrink factor')), )) RPROF = OrderedDict(( ('plot', Conf(None, True, 'o', {'nargs': '?', 'const': '', 'type': str}, False, 'specify which variable to plot')), ('plot_grid', Conf(True, False, None, {}, True, 'plot grid')), ('plot_temperature', Conf(True, False, None, {}, True, 'plot temperature')), ('plot_minmaxtemp', Conf(False, False, None, {}, True, 'plot min and max temperature')), ('plot_velocity', Conf(True, False, None, {}, True, 'plot velocity')), ('plot_minmaxvelo', Conf(False, False, None, {}, True, 'plot min and max velocity')), ('plot_viscosity', Conf(False, False, None, {}, True, 'plot viscosity')), ('plot_minmaxvisco', Conf(False, False, None, {}, True, 'plot min and max viscosity')), ('plot_advection', Conf(True, False, None, {}, True, 'plot heat advction')), ('plot_energy', Conf(True, False, None, {}, True, 'plot energy')), ('plot_concentration', Conf(True, False, None, {}, True, 'plot concentration')), ('plot_minmaxcon', Conf(False, False, None, {}, True, 'plot min and max concentration')), ('plot_conctheo', Conf(True, False, None, {}, True, 'plot concentration theo')), ('plot_overturn_init', Conf(True, False, None, {}, True, 'plot overturn init')), ('plot_difference', Conf(True, False, None, {}, True, 'plot difference between T and C profs and overturned \ version of their initial values')), )) TIME = OrderedDict(( ('compstat', Conf(True, True, None, {}, True, 'compute steady state statistics')), ('annottmin', Conf(False, True, None, {}, True, 'put an arrow at tminc and tmint')), ('tmint', Conf(0., True, None, {}, False, 'specify tmint')), ('tminc', Conf(0., True, None, {}, False, 'specify tminc')), )) PLATES = OrderedDict(( ('vzcheck', Conf(False, True, None, {}, True, 'activate Colin\'s version with vz checking')), ('timeprofile', Conf(False, True, None, {}, True, 'plots nb of plates in function of time')), ('dsa', Conf(0.05, False, None, {}, True, 'thickness of sticky air')), )) VAR = OrderedDict(( )) CONFIG = OrderedDict(( ('create', Conf(None, True, None, {'action': 'store_true'}, False, 'create new config file')), ('update', Conf(None, True, None, {'action': 'store_true'}, False, 'add missing entries to existing config file')), ('edit', Conf(None, True, None, {'action': 'store_true'}, False, 'open config file in a text editor')), ('editor', Conf('vim', False, None, {}, True, 'text editor')), )) def config_cmd(args): """sub command config""" if args.create or args.update: create_config() if args.edit: call(shlex.split(args.editor + ' ' + CONFIG_FILE)) Sub = namedtuple('Sub', ['conf_dict', 'use_core', 'func', 'help_string']) SUB_CMDS = OrderedDict(( ('field', Sub(FIELD, True, commands.field_cmd, 'plot scalar fields')), ('rprof', Sub(RPROF, True, commands.rprof_cmd, 'plot radial profiles')), ('time', Sub(TIME, True, commands.time_cmd, 'plot temporal series')), ('plates', Sub(PLATES, True, commands.plates_cmd, 'plate analysis')), ('var', Sub(VAR, False, commands.var_cmd, 'print the list of variables')), ('config', Sub(CONFIG, False, config_cmd, 'configuration handling')), )) DummySub = namedtuple('DummySub', ['conf_dict']) DUMMY_CMDS = OrderedDict(( ('core', DummySub(CORE)), )) DUMMY_CMDS.update(SUB_CMDS) def _set_conf_default(conf_dict, opt, dflt): """set default value of option in conf_dict""" conf_dict[opt] = conf_dict[opt]._replace(default=dflt) def _read_section(config_parser, sub_cmd, meta): """read section of config parser read section corresponding to the sub command sub_cmd and set meta.conf_dict default values to the read values """ missing_opts = [] for opt, meta_opt in meta.conf_dict.items(): if not config_parser.has_option(sub_cmd, opt): if meta_opt.conf_arg: missing_opts.append(opt) continue if isinstance(meta_opt.default, bool): dflt = config_parser.getboolean(sub_cmd, opt) elif isinstance(meta_opt.default, float): dflt = config_parser.getfloat(sub_cmd, opt) elif isinstance(meta_opt.default, int): dflt = config_parser.getint(sub_cmd, opt) else: dflt = config_parser.get(sub_cmd, opt) _set_conf_default(meta.conf_dict, opt, dflt) return missing_opts def create_config(): """Create config file""" config_parser = configparser.ConfigParser() for sub_cmd, meta in DUMMY_CMDS.items(): config_parser.add_section(sub_cmd) for opt, opt_meta in meta.conf_dict.items(): if opt_meta.conf_arg: config_parser.set(sub_cmd, opt, str(opt_meta.default)) with open(CONFIG_FILE, 'w') as out_stream: config_parser.write(out_stream) def read_config(args): """Read config file and set conf_dict as needed""" if not os.path.isfile(CONFIG_FILE): if not args.update: print('Config file {} not found.'.format(CONFIG_FILE)) print('Run stagpy config --create') return config_parser = configparser.ConfigParser() config_parser.read(CONFIG_FILE) missing_sections = [] for sub_cmd, meta in DUMMY_CMDS.items(): if not config_parser.has_section(sub_cmd): missing_sections.append(sub_cmd) continue missing_opts = _read_section(config_parser, sub_cmd, meta) if missing_opts and not args.update: print('WARNING! Missing options in {} section of config file:'. format(sub_cmd)) print(*missing_opts) print() if missing_sections and not args.update: print('WARNING! Missing sections in config file:') print(*missing_sections) print() if (missing_sections or missing_opts) and not args.update: print('Run stagpy config --update to update config file') print() class Toggle(argparse.Action): """argparse Action to store True/False to a +/-arg""" def __call__(self, parser, namespace, values, option_string=None): """set args attribute with True/False""" setattr(namespace, self.dest, bool('-+'.index(option_string[0]))) def add_args(parser, conf_dict): """Add arguments to a parser""" for arg, conf in conf_dict.items(): if not conf.cmd_arg: continue if isinstance(conf.default, bool): conf.kwargs.update(action=Toggle, nargs=0) names = ['-{}'.format(arg), '+{}'.format(arg)] if conf.shortname is not None: names.append('-{}'.format(conf.shortname)) names.append('+{}'.format(conf.shortname)) else: if conf.default is not None: conf.kwargs.setdefault('type', type(conf.default)) names = ['--{}'.format(arg)] if conf.shortname is not None: names.append('-{}'.format(conf.shortname)) conf.kwargs.update(help=conf.help_string) parser.add_argument(*names, **conf.kwargs) parser.set_defaults(**{a: c.default for a, c in conf_dict.items()}) return parser def parse_args(): """Parse cmd line arguments""" # get path from config file before if not os.path.isdir(CONFIG_DIR): mkdir(CONFIG_DIR) dummy_parser = argparse.ArgumentParser(add_help=False) _, remainder = dummy_parser.parse_known_args() keep_cmd_path = '-p' in remainder or '--path' in remainder dummy_parser = add_args(dummy_parser, {'path': CORE['path']}) args, remainder = dummy_parser.parse_known_args() cmd_path = args.path _set_conf_default(CORE, 'path', args.path) if 'config' in remainder: dummy_sub = dummy_parser.add_subparsers() dummy_conf = dummy_sub.add_parser('config', add_help=False) dummy_conf = add_args(dummy_conf, CONFIG) args, _ = dummy_parser.parse_known_args() else: args.create = False args.edit = False args.update = False if not (args.create or args.edit): try: read_config(args) except: print('ERROR while reading config file') print('Run stagpy config --create to obtain a new config file') print('=' * 26) raise if keep_cmd_path: args.path = cmd_path _set_conf_default(CORE, 'path', args.path) par_nml = parfile.readpar(args) main_parser = argparse.ArgumentParser( description='read and process StagYY binary data') main_parser = add_args(main_parser, {'path': CORE['path']}) main_parser.set_defaults(func=lambda _: print('stagpy -h for usage')) subparsers = main_parser.add_subparsers() core_parser = argparse.ArgumentParser(add_help=False, prefix_chars='-+') core_parser = add_args(core_parser, CORE) core_parser.set_defaults(name=par_nml['ioin']['output_file_stem']) for sub_cmd, meta in SUB_CMDS.items(): kwargs = {'prefix_chars': '+-', 'help': meta.help_string} if meta.use_core: kwargs.update(parents=[core_parser]) dummy_parser = subparsers.add_parser(sub_cmd, **kwargs) dummy_parser = add_args(dummy_parser, meta.conf_dict) dummy_parser.set_defaults(func=meta.func) argcomplete.autocomplete(main_parser) args = main_parser.parse_args() args.par_nml = par_nml return args PKXILH=G~ stagpy/field.py"""plot fields""" import numpy as np from . import constants, misc from .stagdata import BinData def plot_scalar(args, stgdat, var): """var: one of the key of constants.FIELD_VAR_LIST""" plt = args.plt if var == 's': fld = stgdat.calc_stream() else: fld = stgdat.fields[var] # adding a row at the end to have continuous field if stgdat.geom == 'annulus': if stgdat.par_type == 'vp': if var != 's': fld = fld[:, :, 0].T else: newline = fld[:, 0, 0] fld = np.vstack([fld[:, :, 0].T, newline]) xmesh, ymesh = stgdat.x_mesh[0, :, :], stgdat.y_mesh[0, :, :] fig, axis = plt.subplots(ncols=1) if stgdat.geom == 'annulus': if var == 'n': surf = axis.pcolormesh(xmesh, ymesh, fld, norm=args.mpl.colors.LogNorm(), cmap='jet_r', rasterized=not args.pdf, shading='gouraud') elif var == 'd': surf = axis.pcolormesh(xmesh, ymesh, fld, cmap='bwr_r', vmin=0.96, vmax=1.04, rasterized=not args.pdf, shading='gouraud') else: surf = axis.pcolormesh(xmesh, ymesh, fld, cmap='jet', rasterized=not args.pdf, shading='gouraud') cbar = plt.colorbar(surf, shrink=args.shrinkcb) cbar.set_label(constants.FIELD_VAR_LIST[var].name) plt.axis('equal') plt.axis('off') return fig, axis def plot_stream(args, fig, axis, component1, component2): """use of streamplot to plot stream lines only works in cartesian with regular grids """ x_1, v_1 = component1 x_2, v_2 = component2 v_tot = np.sqrt(v_1**2 + v_2**2) lwd = 2 * v_tot / v_tot.max() args.plt.figure(fig.number) axis.streamplot(x_1, x_2, v_1, v_2, density=0.8, color='k', linewidth=lwd) def field_cmd(args): """extract and plot field data""" for timestep in range(*args.timestep): for var, meta in constants.FIELD_VAR_LIST.items(): if misc.get_arg(args, meta.arg): # will read vp many times! stgdat = BinData(args, var, timestep) fig, _ = plot_scalar(args, stgdat, var) args.plt.figure(fig.number) args.plt.tight_layout() args.plt.savefig( misc.out_name(args, var).format(stgdat.step) + '.pdf', format='PDF') args.plt.close(fig) PKXILHqqhQ''stagpy/stagdata.py"""define StagyyData""" import numpy as np import re import struct from itertools import zip_longest from scipy import integrate from . import constants, misc class BinData: """reads StagYY binary data and processes them""" def __init__(self, args, var, timestep): """read the necessary binary file after init, the StagyyData object is ready for processing """ self.args = args self.var = var self.par_type = constants.FIELD_VAR_LIST[var].par self.geom = args.geometry self.file_format = 'l' self.step = timestep # name of the file to read self.fullname = misc.stag_file(args, self.par_type, timestep) self.nval = 4 if self.par_type == 'vp' else 1 with open(self.fullname, 'rb') as self._fid: self._catch_header() self._readfile() def _readbin(self, fmt='i', nwords=1, nbytes=4): """read n words of n bytes with fmt format Return a tuple of elements if more than one element. Default: read 1 word of 4 bytes formatted as an integer. """ elts = struct.unpack(fmt * nwords, self._fid.read(nwords * nbytes)) if len(elts) == 1: elts = elts[0] return elts def _catch_header(self): """reads header of binary file""" self.nmagic = self._readbin() # Version # check nb components if (self.nmagic < 100 and self.nval > 1) \ or (self.nmagic > 300 and self.nval == 1): raise ValueError('wrong number of components in field') # extra ghost point in horizontal direction self.xyp = int((self.nmagic % 100) >= 9 and self.nval == 4) # total number of values in # latitude, longitude and radius directions self.nthtot, self.nphtot, self.nrtot = self._readbin(nwords=3) # number of blocks, 2 for yinyang self.nblocks = self._readbin() # Aspect ratio self.aspect = self._readbin('f', 2) self.aspect = np.array(self.aspect) # Number of parallel subdomains in the th,ph,r and b directions self.nnth, self.nnph, self.nnr = self._readbin(nwords=3) self.nnb = self._readbin() self.nr2 = self.nrtot * 2 + 1 self.rgeom = self._readbin('f', self.nr2) # r-coordinates self.rcmb = self._readbin('f') # radius of the cmb self.ti_step = self._readbin() self.ti_ad = self._readbin('f') self.erupta_total = self._readbin('f') self.bot_temp = self._readbin('f') # theta coordinates self.th_coord = np.array(self._readbin('f', self.nthtot)) # force to pi/2 if 2D self.th_coord = np.array(np.pi / 2) # phi coordinates ph_coord = np.array(self._readbin('f', self.nphtot)) self._ph_coord = ph_coord # to have continuous field self.ph_coord = np.append(ph_coord, ph_coord[1] - ph_coord[0]) # radius coordinates self.r_coord = np.array(self._readbin('f', self.nrtot)) # create meshgrids self.th_mesh, self.ph_mesh, self.r_mesh = np.meshgrid( self.th_coord, self.ph_coord, self.r_coord + self.rcmb, indexing='ij') # compute cartesian coordinates # z along rotation axis at theta=0 # x at th=90, phi=0 # y at th=90, phi=90 self.x_mesh = self.r_mesh * np.cos(self.ph_mesh) * np.sin(self.th_mesh) self.y_mesh = self.r_mesh * np.sin(self.ph_mesh) * np.sin(self.th_mesh) self.z_mesh = self.r_mesh * np.cos(self.th_mesh) def _readfile(self): """read scalar/vector fields""" # compute nth, nph, nr and nb PER CPU nth = self.nthtot // self.nnth nph = self.nphtot // self.nnph nrd = self.nrtot // self.nnr nbk = self.nblocks // self.nnb # the number of values per 'read' block npi = (nth + self.xyp) * (nph + self.xyp) * nrd * nbk * self.nval if self.nval > 1: self.scalefac = self._readbin('f') else: self.scalefac = 1 # flds should be construct with the "normal" indexing order th, ph, r # there shouldn't be a need to transpose in plot_scalar dim_fields = (self.nblocks, self.nrtot, self.nphtot + self.xyp, self.nthtot + self.xyp) flds = [] for _ in range(self.nval): flds.append(np.zeros(dim_fields)) # loop over parallel subdomains for ibc in np.arange(self.nnb): for irc in np.arange(self.nnr): for iphc in np.arange(self.nnph): for ithc in np.arange(self.nnth): # read the data for this CPU file_content = self._readbin('f', npi) data_cpu = np.array(file_content) * self.scalefac # Create a 3D matrix from these data data_cpu_3d = data_cpu.reshape( (nbk, nrd, nph + self.xyp, nth + self.xyp, self.nval)) # Add local 3D matrix to global matrix sth = ithc * nth eth = ithc * nth + nth + self.xyp sph = iphc * nph eph = iphc * nph + nph + self.xyp srd = irc * nrd erd = irc * nrd + nrd snb = ibc * nbk enb = ibc * nbk + nbk for idx, fld in enumerate(flds): fld[snb:enb, srd:erd, sph:eph, sth:eth] = \ data_cpu_3d[:, :, :, :, idx] self.fields = {} fld_names = ['u', 'v', 'w', 'p'] if self.par_type == 'vp' \ else [self.var] for fld_name, fld in zip(fld_names, flds): self.fields[fld_name] = fld[0, :, :, :] def calc_stream(self): """computes and returns the stream function only make sense with vp fields """ # should add test if vp fields or not vphi = self.fields['v'][:, :, 0] # interpolate to the same phi vph2 = -0.5 * (vphi + np.roll(vphi, 1, 1)) v_r = self.fields['w'][:, :, 0] n_r, nph = np.shape(v_r) stream = np.zeros(np.shape(vphi)) # integrate first on phi stream[0, 1:nph - 1] = self.rcmb * \ integrate.cumtrapz(v_r[0, 0:nph - 1], self._ph_coord) stream[0, 0] = 0 # use r coordinates where vphi is defined rcoord = self.rcmb + np.array( self.rgeom[0:np.shape(self.rgeom)[0] - 1:2]) for iph in range(0, np.shape(vph2)[1] - 1): stream[1:n_r, iph] = stream[0, iph] + \ integrate.cumtrapz(vph2[:, iph], rcoord) # integrate on r stream = stream - np.mean(stream[n_r / 2, :]) # remove some typical value. # Would be better to compute the global average # taking into account variable grid spacing return stream class RprofData: """extract radial profiles data""" def __init__(self, args): """create RprofData object""" step_regex = re.compile(r'^\*+step:\s*(\d+) ; time =\s*(\S+)') self._readproffile(args, step_regex) def _readproffile(self, args, step_regex): """extract info from rprof.dat""" proffile = misc.stag_file(args, 'rprof.dat') timesteps = [] data0 = [] lnum = -1 with open(proffile) as stream: for line in stream: if line != '\n': lnum += 1 if line[0] == '*': match = step_regex.match(line) timesteps.append([lnum, int(match.group(1)), float(match.group(2))]) else: data0.append(np.array(line.split())) tsteps = np.array(timesteps) nsteps = tsteps.shape[0] data = np.array(data0) # all the processings of timesteps # should be in commands.*_cmd # instead of main.py # since it could be different between # the different modules istart, ilast, istep = args.timestep if ilast == -1: ilast = nsteps - 1 if istart == -1: istart = nsteps - 1 args.timestep = istart, ilast, istep nzp = [] for iti in range(0, nsteps - 1): nzp = np.append(nzp, tsteps[iti + 1, 0] - tsteps[iti, 0] - 1) nzp = np.append(nzp, lnum - tsteps[nsteps - 1, 0]) nzs = [[0, 0, 0]] nzc = 0 for iti in range(1, nsteps): if nzp[iti] != nzp[iti - 1]: nzs.append([iti, iti - nzc, int(nzp[iti - 1])]) nzc = iti if nzp[nsteps - 1] != nzs[-1][1]: nzs.append([nsteps, nsteps - nzc, int(nzp[nsteps - 1])]) nzi = np.array(nzs) self.data = data self.tsteps = tsteps self.nzi = nzi class TimeData: """extract temporal series""" def __init__(self, args): """read temporal series from time.dat""" timefile = misc.stag_file(args, 'time.dat') with open(timefile, 'r') as infile: first = infile.readline() line = infile.readline() data = [list(misc.parse_line(line))] for line in infile: step = list(misc.parse_line(line, convert=[int])) # remove useless lines when run is restarted while step[0] <= data[-1][0]: data.pop() data.append(step) self.colnames = first.split() # suppress two columns from the header. # Only temporary since this has been corrected in stag # WARNING: possibly a problem is some columns are added? if len(self.colnames) == 33: self.colnames = self.colnames[:28] + self.colnames[30:] self.data = np.array(list(zip_longest(*data, fillvalue=0))).T PKI3Hstagpy/__init__.pyPKXILHoq5I5Istagpy/plates.py"""Plots position of subduction and ridge at the surface. Date: 2016/26/01 """ import numpy as np import sys from . import misc from .stagdata import BinData, RprofData, TimeData from scipy.signal import argrelextrema from copy import deepcopy def detect_plates_vzcheck(stagdat_t, stagdat_vp, stagdat_h, rprof_data, args, seuil_memz): """detect plates and check with vz and plate size""" v_z = stagdat_vp.fields['w'] v_x = stagdat_vp.fields['v'] h2o = stagdat_h.fields['h'] tcell = stagdat_t.fields['t'] data = rprof_data.data n_z = len(v_z) nphi = len(v_z[0]) - 1 radius = list(map(float, data[0:n_z, 0])) if args.par_nml['geometry']['shape'].lower() == 'spherical': rcmb = args.par_nml['geometry']['r_cmb'] else: rcmb = 0. dphi = 1 / nphi # calculing radius on the grid radiusgrid = len(radius) * [0] radiusgrid.append(1) for i in range(1, len(radius)): radiusgrid[i] = 2 * radius[i - 1] - radiusgrid[i - 1] for i in range(len(radiusgrid)): radiusgrid[i] += rcmb for i in range(len(radius)): radius[i] += rcmb # water profile water_profile = n_z * [0] for i_z in range(n_z): for phi in range(nphi): water_profile[i_z] += h2o[i_z, phi, 0] / nphi # calculing tmean tmean = 0 for i_r in range(len(radius)): for phi in range(nphi): tmean += (radiusgrid[i_r + 1]**2 - radiusgrid[i_r] ** 2) * dphi * tcell[i_r, phi] tmean /= (radiusgrid[-1]**2 - rcmb**2) # calculing temperature on the grid and vz_mean/v_rms v_rms = 0 vz_mean = 0 tgrid = np.zeros((n_z + 1, nphi)) for phi in range(nphi): tgrid[0, phi] = 1 for i_z in range(1, n_z): for phi in range(nphi): tgrid[i_z, phi] = ( tcell[i_z - 1, phi] * (radiusgrid[i_z] - radius[i_z - 1]) + tcell[i_z, phi] * (-radiusgrid[i_z] + radius[i_z])) / (radius[i_z] - radius[i_z - 1]) v_rms += (v_z[i_z, phi, 0]**2 + v_x[i_z, phi, 0]**2) / (nphi * n_z) vz_mean += abs(v_z[i_z, phi, 0]) / (nphi * n_z) v_rms = v_rms**0.5 print(v_rms, vz_mean) flux_c = n_z * [0] for i_z in range(1, n_z - 1): for phi in range(nphi): flux_c[i_z] += (tgrid[i_z, phi] - tmean) * \ v_z[i_z, phi, 0] * radiusgrid[i_z] * dphi # checking stagnant lid stagnant_lid = True max_flx = np.max(flux_c) for i_z in range(n_z - n_z // 20, n_z): if abs(flux_c[i_z]) > max_flx / 50: stagnant_lid = False break if stagnant_lid: print('stagnant lid') sys.exit() else: # verifying horizontal plate speed and closeness of plates dvphi = nphi * [0] dvx_thres = 16 * v_rms for phi in range(0, nphi): dvphi[phi] = (v_x[n_z - 1, phi, 0] - v_x[n_z - 1, phi - 1, 0]) / ((1 + rcmb) * dphi) limits = [] for phi in range(0, nphi - nphi // 33): mark = True for i in range(phi - nphi // 33, phi + nphi // 33): if abs(dvphi[i]) > abs(dvphi[phi]): mark = False if mark and abs(dvphi[phi]) >= dvx_thres: limits.append(phi) for phi in range(nphi - nphi // 33 + 1, nphi): mark = True for i in range(phi - nphi // 33 - nphi, phi + nphi // 33 - nphi): if abs(dvphi[i]) > abs(dvphi[phi]): mark = False if mark and abs(dvphi[phi]) >= dvx_thres: limits.append(phi) print(limits) # verifying vertical speed k = 0 for i in range(len(limits)): vzm = 0 phi = limits[i - k] if phi == nphi - 1: for i_z in range(1, n_z): vzm += (abs(v_z[i_z, phi, 0]) + abs(v_z[i_z, phi - 1, 0]) + abs(v_z[i_z, 0, 0])) / (n_z * 3) else: for i_z in range(0, n_z): vzm += (abs(v_z[i_z, phi, 0]) + abs(v_z[i_z, phi - 1, 0]) + abs(v_z[i_z, phi + 1, 0])) / (n_z * 3) if seuil_memz != 0: vz_thres = vz_mean * 0.1 + seuil_memz / 2 else: vz_thres = vz_mean * 0 if vzm < vz_thres: limits.remove(phi) k += 1 print(limits) print('\n') return limits, nphi, dvphi, vz_thres, v_x[n_z - 1, :, 0], water_profile def detect_plates(args, velocity): """detect plates using horizontal velocity""" velocityfld = velocity.fields['v'] ph_coord = velocity.ph_coord dsa = args.dsa # we are a bit below the surface; should check if you are in the # mechanical/thermal boundary layer indsurf = np.argmin(abs((1 - dsa) - velocity.r_coord)) - 4 vphi = velocityfld[:, :, 0] vph2 = 0.5 * (vphi + np.roll(vphi, 1, 1)) # interpolate to the same phi # velocity derivation dvph2 = (np.diff(vph2[indsurf, :]) / (ph_coord[0] * 2.)) # prepare stuff to find trenches and ridges myorder_trench = 40 myorder_ridge = 20 # threshold # finding trenches pom2 = deepcopy(dvph2) maskbigdvel = np.amin(dvph2) * 0.25 # putting threshold pom2[pom2 > maskbigdvel] = 0 # user putting threshold argless_dv = argrelextrema( pom2, np.less, order=myorder_trench, mode='wrap')[0] trench = ph_coord[argless_dv] # finding ridges pom2 = deepcopy(dvph2) masksmalldvel = np.amax(dvph2) * 0.2 # putting threshold pom2[pom2 < masksmalldvel] = 0 arggreat_dv = argrelextrema( pom2, np.greater, order=myorder_ridge, mode='wrap')[0] ridge = ph_coord[arggreat_dv] # agetrench=age_surface_dim[argless_dv] # age of the trench # elimination of ridges that are too close to trench argdel = [] if trench and ridge: for i in range(len(ridge)): mdistance = np.amin(abs(trench - ridge[i])) if mdistance < 0.016: argdel.append(i) if argdel: print('deleting from ridge', trench, ridge[argdel]) ridge = np.delete(ridge, np.array(argdel)) arggreat_dv = np.delete(arggreat_dv, np.array(argdel)) return trench, ridge def plot_plates(args, velocity, temp, conc, age, timestep, trench, ridge): """handle ploting stuffs""" plt = args.plt lwd = args.linewidth meanvrms = 605.0 # to be changed ttransit = 1.78e15 # My yearins = 2.16E7 dsa = 0.05 plot_age = True velocityfld = velocity.fields['v'] tempfld = temp.fields['t'] concfld = conc.fields['c'] agefld = age.fields['a'] # if stgdat.par_type == 'vp': # fld = fld[:, :, 0] newline = tempfld[:, 0, 0] tempfld = np.vstack([tempfld[:, :, 0].T, newline]).T newline = concfld[:, 0, 0] concfld = np.vstack([concfld[:, :, 0].T, newline]).T newline = agefld[:, 0, 0] agefld = np.vstack([agefld[:, :, 0].T, newline]).T # we are a bit below the surface; delete "-some number" to be just below # the surface (that is considered plane here); should check if you are in # the mechanical/thermal boundary layer indsurf = np.argmin(abs((1 - dsa) - temp.r_coord)) - 4 # depth to detect the continents indcont = np.argmin(abs((1 - dsa) - np.array(velocity.r_coord))) - 10 continents = np.ma.masked_where( np.logical_or(concfld[indcont, :-1] < 3, concfld[indcont, :-1] > 4), concfld[indcont, :-1]) # masked array, only continents are true continentsall = continents / continents # if(vp.r_coord[indsurf]>1.-dsa): # print 'WARNING lowering index for surface' # indsurf=indsurf-1 # if verbose_figures: # age just below the surface if plot_age: age_surface = np.ma.masked_where( agefld[indsurf, :] < 0.00001, agefld[indsurf, :]) age_surface_dim = age_surface * meanvrms * ttransit / yearins / 1.e6 ph_coord = conc.ph_coord # velocity vphi = velocityfld[:, :, 0] vph2 = 0.5 * (vphi + np.roll(vphi, 1, 1)) # interpolate to the same phi dvph2 = (np.diff(vph2[indsurf, :]) / (ph_coord[0] * 2.)) # dvph2=dvph2/amax(abs(dvph2)) # normalization # plotting _, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, sharex=True, figsize=(10, 12)) ax1.plot(ph_coord[:-1], concfld[indsurf, :-1], color='g', linewidth=lwd, label='Conc') ax2.plot(ph_coord[:-1], tempfld[indsurf, :-1], color='m', linewidth=lwd, label='Temp') ax3.plot(ph_coord[:-1] + ph_coord[0], dvph2, color='c', linewidth=lwd, label='dv') ax4.plot(ph_coord[:-1], vph2[indsurf, :-1], linewidth=lwd, label='Vel') ax1.fill_between( ph_coord[:-1], continents, 1., facecolor='gray', alpha=0.25) ax2.fill_between( ph_coord[:-1], continentsall, 0., facecolor='gray', alpha=0.25) ax2.set_ylim(0, 1) ax3.fill_between( ph_coord[:-1], continentsall * round(1.5 * np.amax(dvph2), 1), round(np.amin(dvph2) * 1.1, 1), facecolor='gray', alpha=0.25) ax3.set_ylim( round(np.amin(dvph2) * 1.1, 1), round(1.5 * np.amax(dvph2), 1)) ax4.fill_between( ph_coord[:-1], continentsall * 5e3, -5000, facecolor='gray', alpha=0.25) ax4.set_ylim(-5000, 5000) ax1.set_ylabel("Concentration") ax2.set_ylabel("Temperature") ax3.set_ylabel("dv") ax4.set_ylabel("Velocity") ax1.set_title(timestep) # topography fname = misc.stag_file(args, 'sc', timestep=temp.step, suffix='.dat') depth_mantle = 2890.0 # in km topo = np.genfromtxt(fname) # rescaling topography! topo[:, 1] = topo[:, 1] / (1. - dsa) topomin = -50 topomax = 50 # majorLocator = MultipleLocator(20) ax31 = ax3.twinx() ax31.set_ylabel("Topography [km]") ax31.plot(topo[:, 0], topo[:, 1] * depth_mantle, color='black', alpha=0.4) ax31.set_ylim(topomin, topomax) ax41 = ax4.twinx() ax41.set_ylabel("Topography [km]") ax41.axhline( y=0, xmin=0, xmax=2 * np.pi, color='black', ls='dashed', alpha=0.7) for i in range(len(trench)): ax41.axvline( x=trench[i], ymin=topomin, ymax=topomax, color='red', ls='dashed', alpha=0.4) for i in range(len(ridge)): ax41.axvline( x=ridge[i], ymin=topomin, ymax=topomax, color='green', ls='dashed', alpha=0.8) ax41.plot(topo[:, 0], topo[:, 1] * depth_mantle, color='black', alpha=0.7) ax41.set_ylim(topomin, topomax) ax1.set_xlim(0, 2 * np.pi) figname = misc.out_name(args, 'surf').format(temp.step) + '.pdf' plt.savefig(figname, format='PDF') plt.close() # plotting only velocity and topography _, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(12, 8)) ax1.plot(ph_coord[:-1], vph2[indsurf, :-1], linewidth=lwd, label='Vel') ax1.axhline(y=0, xmin=0, xmax=2 * np.pi, color='black', ls='dashed', alpha=0.4) ax1.set_ylim(-5000, 5000) ax1.set_ylabel("Velocity") topomax = 30 topomin = -60 for i in range(len(trench)): ax1.axvline( x=trench[i], ymin=topomin, ymax=topomax, color='red', ls='dashed', alpha=0.8) # detection of the distance in between subduction and continent ph_coord_noendpoint = ph_coord[:-1] distancecont = min( abs(ph_coord_noendpoint[continentsall == 1] - trench[i])) argdistancecont = np.argmin( abs(ph_coord_noendpoint[continentsall == 1] - trench[i])) continentpos = ph_coord_noendpoint[continentsall == 1][argdistancecont] # do i have a ridge in between continent edge and trench? if ridge: if min(abs(continentpos - ridge)) > distancecont: # unexistent variables? # ph_trench_subd.append(trench[i]) # age_subd.append(agetrench[i]) # ph_cont_subd.append(continentpos) # distance_subd.append(distancecont) # times_subd.append(temp.ti_ad) # continent is on the left if (continentpos - trench[i]) < 0: ax1.annotate('', xy=(trench[i] - distancecont, 2000), xycoords='data', xytext=(trench[i], 2000), textcoords='data', arrowprops=dict(arrowstyle="->", shrinkA=0, shrinkB=0)) else: # continent is on the right ax1.annotate('', xy=(trench[i] + distancecont, 2000), xycoords='data', xytext=(trench[i], 2000), textcoords='data', arrowprops=dict(arrowstyle="->", shrinkA=0, shrinkB=0)) ax1.axvline( x=trench[i], ymin=topomin, ymax=topomax, color='red', ls='dashed', alpha=0.8) ax1.grid() for i in range(len(ridge)): ax1.axvline( x=ridge[i], ymin=topomin, ymax=topomax, color='green', ls='dashed', alpha=0.8) ax2.set_ylabel("Topography [km]") ax2.plot(topo[:, 0], topo[:, 1] * depth_mantle, color='black') ax2.set_xlim(0, 2 * np.pi) dtopo = deepcopy(topo[:, 1] * depth_mantle) mask = dtopo > 0 water = deepcopy(dtopo) water[mask] = 0 ax2.set_ylim(topomin, topomax) for i in range(len(trench)): ax2.axvline( x=trench[i], ymin=topomin, ymax=topomax, color='red', ls='dashed', alpha=0.8) for i in range(len(ridge)): ax2.axvline( x=ridge[i], ymin=topomin, ymax=topomax, color='green', ls='dashed', alpha=0.8) ax1.set_title(timestep) figname = misc.out_name(args, 'surfvel').format(temp.step) + '.pdf' plt.savefig(figname, format='PDF') plt.close() # plotting only velocity and age at surface if plot_age: _, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(12, 8)) ax1.plot(ph_coord[:-1], vph2[indsurf, :-1], linewidth=lwd, label='Vel') ax1.axhline( y=0, xmin=0, xmax=2 * np.pi, color='black', ls='dashed', alpha=0.4) ax1.set_ylim(-5000, 5000) ax1.set_ylabel("Velocity") agemax = 280 agemin = 0 for i in range(len(trench)): ax1.axvline( x=trench[i], ymin=agemin, ymax=agemax, color='red', ls='dashed', alpha=0.8) for i in range(len(ridge)): ax1.axvline( x=ridge[i], ymin=agemin, ymax=agemax, color='green', ls='dashed', alpha=0.8) ax2.set_ylabel("Age [My]") # in dimensions ax2.plot(ph_coord[:-1], age_surface_dim[:-1], color='black') ax2.set_xlim(0, 2 * np.pi) ax2.fill_between( ph_coord[:-1], continentsall * agemax, agemin, facecolor='#8B6914', alpha=0.2) ax2.set_ylim(agemin, agemax) for i in range(len(trench)): ax2.axvline( x=trench[i], ymin=agemin, ymax=agemax, color='red', ls='dashed', alpha=0.8) for i in range(len(ridge)): ax2.axvline( x=ridge[i], ymin=agemin, ymax=agemax, color='green', ls='dashed', alpha=0.8) ax1.set_title(timestep) figname = misc.out_name(args, 'surfage').format(temp.step) + '.pdf' plt.savefig(figname, format='PDF') plt.close() return None def plates_cmd(args): """find positions of trenches and subductions uses velocity field (velocity derivation) plots the number of plates over a designated lapse of time """ if args.vzcheck: seuil_memz = 0 nb_plates = [] timedat = TimeData(args) slc = slice(*(i * args.par_nml['ioin']['save_file_framestep'] for i in args.timestep)) time, ch2o = timedat.data[:, 1][slc], timedat.data[:, 27][slc] for timestep in range(*args.timestep): velocity = BinData(args, 'v', timestep) temp = BinData(args, 't', timestep) if args.vzcheck: rprof_data = RprofData(args) water = BinData(args, 'h', timestep) rprof_data = RprofData(args) plt = args.plt limits, nphi, dvphi, seuil_memz, vphi_surf, water_profile =\ detect_plates_vzcheck(temp, velocity, water, rprof_data, args, seuil_memz) limits.sort() sizeplates = [limits[0] + nphi - limits[-1]] for lim in range(1, len(limits)): sizeplates.append(limits[lim] - limits[lim - 1]) lim = len(limits) * [max(dvphi)] plt.figure(timestep) plt.subplot(221) plt.axis([0, len(velocity.fields['w'][0]) - 1, np.min(vphi_surf) * 1.2, np.max(vphi_surf) * 1.2]) plt.plot(vphi_surf) plt.subplot(223) plt.axis( [0, len(velocity.fields['w'][0]) - 1, np.min(dvphi) * 1.2, np.max(dvphi) * 1.2]) plt.plot(dvphi) plt.scatter(limits, lim, color='red') plt.subplot(222) plt.hist(sizeplates, 10, (0, nphi / 2)) plt.subplot(224) plt.plot(water_profile) plt.savefig('plates' + str(timestep) + '.pdf', format='PDF') nb_plates.append(len(limits)) plt.close(timestep) else: conc = BinData(args, 'c', timestep) age = BinData(args, 'a', timestep) trenches, ridges = detect_plates(args, velocity) plot_plates( args, velocity, temp, conc, age, timestep, trenches, ridges) if args.timeprofile and args.vzcheck: for i in range(2, len(nb_plates) - 3): nb_plates[i] = (nb_plates[i - 2] + nb_plates[i - 1] + nb_plates[i] + nb_plates[i + 1] + nb_plates[i + 2]) / 5 plt.figure(-1) plt.subplot(121) plt.axis([time[0], time[-1], 0, np.max(nb_plates)]) plt.plot(time, nb_plates) plt.subplot(122) plt.plot(time, ch2o) plt.savefig('plates_{}_{}_{}.pdf'.format(*args.timestep), format='PDF') plt.close(-1) PKXILHGc;;stagpy/rprof.py"""Plots radial profiles coming out of stagyy. Author: Stephane Labrosse with inputs from Martina Ulvrova and Adrien Morison Date: 2015/09/11 """ import numpy as np from scipy import integrate as itg import math from . import constants, misc from .stagdata import RprofData def _normprof(rrr, func): # for args.plot_difference """Volumetric norm of a profile Two arrays: rrr is the radius position and f the function. """ norm = 3. / (rrr[-1]**3 - rrr[0]**3) * itg.trapz(func**2 * rrr**2, rrr) return norm def _extrap(xpos, xpoints, ypoints): # for args.plot_difference """np.interp function with linear extrapolation. Would be best to use degree 3 extrapolation """ ypos = np.interp(xpos, xpoints, ypoints) ypos[xpos < xpoints[0]] = ypoints[0]\ + (xpos[xpos < xpoints[0]] - xpoints[0])\ * (ypoints[0] - ypoints[1]) / (xpoints[0] - xpoints[1]) ypos[xpos > xpoints[-1]] = ypoints[-1]\ + (xpos[xpos > xpoints[-1]] - xpoints[-1])\ * (ypoints[-1] - ypoints[-2]) / (xpoints[-1] - xpoints[-2]) return ypos def _calc_energy(data, ir0, ir1): # for args.plot_energy """Compute energy balance(r)""" zgrid = np.array(data[ir0:ir1, 63], float) zgrid = np.append(zgrid, 1.) dzg = np.array(data[ir0 + 1:ir1, 0], float)\ - np.array(data[ir0:ir1 - 1, 0], float) qadv = np.array(data[ir0:ir1 - 1, 60], float) qadv = np.insert(qadv, 0, 0.) qadv = np.append(qadv, 0.) qcond = (np.array(data[ir0:ir1 - 1, 1], float) - np.array(data[ir0 + 1:ir1, 1], float)) / dzg qcond0 = (1. - float(data[ir0, 1])) / float(data[ir0, 0]) qtop = float(data[ir1, 1]) / (1. - float(data[ir1, 0])) qcond = np.insert(qcond, 0, qcond0) qcond = np.append(qcond, qtop) qtot = qadv + qcond return qtot, qadv, qcond, zgrid def plotprofiles(quant, vartuple, data, tsteps, nzi, rbounds, args, ctheoarg, integrate=False): """Plot the chosen profiles for the chosen timesteps quant holds the strings for the x axis annotation and the legends for the additional profiles vartuple contains the numbers of the column to be plotted """ plt = args.plt istart, ilast, istep = args.timestep lwdth = args.linewidth ftsz = args.fontsize rmin, rmax, rcmb = rbounds axax, initprof = ctheoarg linestyles = ('-', '--', '-.', ':') if integrate: def integ(fct, rad): """(theta, phi) surface scaling factor""" return fct * (rad / rmax)**2 if quant[0] == 'Grid': fig, axe = plt.subplots(2, sharex=True) else: fig, axe = plt.subplots() timename = str(istart) + "_" + str(ilast) + "_" + str(istep) if args.plot_difference: concdif = [] tempdif = [] wmax = [] for step in range(istart + 1, ilast + 1, istep): # find the indices ann = sorted(np.append(nzi[:, 0], step)) inn = ann.index(step) nnz = np.multiply(nzi[:, 1], nzi[:, 2]) ir0 = np.sum([nnz[0:inn]]) + (step - nzi[inn - 1, 0] - 1) * nzi[inn, 2] ir1 = ir0 + nzi[inn, 2] - 1 if quant[0] == 'Energy': energy = _calc_energy(data, ir0, ir1) # Plot the profiles if quant[0] == 'Grid': axe[0].plot(data[ir0:ir1, 0], '-ko', label='z') axe[0].set_ylabel('z', fontsize=ftsz) dzgrid = (np.array(data[ir0 + 1:ir1, 0], np.float) - np.array(data[ir0:ir1 - 1, 0], np.float)) axe[1].plot(dzgrid, '-ko', label='dz') axe[1].set_xlabel('cell number', fontsize=ftsz) axe[1].set_ylabel('dz', fontsize=ftsz) else: if quant[0] == 'Energy': profiles = np.array(np.transpose(energy)[:, [0, 1, 2]], float) radius = np.array(np.transpose(energy)[:, 3], float) + rcmb else: profiles = np.array(data[ir0:ir1, vartuple], float) radius = np.array(data[ir0:ir1, 0], float) + rcmb for i in range(profiles.shape[1]): if integrate: donnee = list(map(integ, profiles[:, i], radius)) else: donnee = profiles[:, i] if i == 0: pplot = plt.plot(donnee, radius, linewidth=lwdth, label=r'$t=%.2e$' % (tsteps[step - 1, 2])) # get color and size characteristics col = pplot[0].get_color() # overturned version of the initial profiles if quant[0] in ('Concentration', 'Temperature') and\ (args.plot_overturn_init or args.plot_difference) and\ step == istart + 1: rfin = (rmax**3 + rmin**3 - radius**3)**(1 / 3) if quant[0] == 'Concentration': conc0 = _extrap(rfin, radius, profiles[:, 0]) if quant[0] == 'Temperature': temp0 = _extrap(rfin, radius, profiles[:, 0]) plt.plot(donnee, rfin, '--', c=col, linewidth=lwdth, label='Overturned') if quant[0] == 'Concentration' and args.plot_difference: concd1 = _normprof(radius, profiles[:, 0] - conc0) concdif.append(concd1) if quant[0] == 'Temperature' and args.plot_difference: tempd1 = _normprof(radius, profiles[:, 0] - temp0) tempdif.append(tempd1) wmax.append(max(np.array(data[ir0:ir1, 7], np.float))) # plot the overturned version of the initial profiles # if ((quant[0] == 'Concentration' or # quant[0] == 'Temperature') and # args.plot_overturn_init and step == istart+1): # rfin = (rmax**3.+rmin**3.-radius**3.)**(1./3.) # plt.plot(donnee, rfin, '--', c=col, # linewidth=lwdth, label='Overturned') # plot the theoretical initial profile and its # overturned version if (quant[0] == 'Concentration' and args.plot_conctheo and step == istart + 1): # plot the full profile between rmin and rmax radius2 = np.linspace(rmin, rmax, 1000) cinit = list(map(initprof, radius2)) rfin = (rmax**3 + rmin**3 - radius2**3)**(1 / 3) plt.plot(cinit, radius2, 'r--', linewidth=lwdth, label='Theoretical') plt.plot(cinit, rfin, 'r-.', linewidth=lwdth, label='Overturned') # add the begining and end points of the stagyy # profile plt.plot([donnee[0], donnee[-1]], [radius[0], radius[-1]], "o", label='StagYY profile ends') plt.xlim([0.9 * donnee[0], 1.2 * donnee[-1]]) else: # additional plots (e. g. min, max) plt.plot(donnee, radius, c=col, dash_capstyle='round', linestyle=linestyles[i], linewidth=lwdth) # change the vertical limits plt.ylim([rmin - 0.05, rmax + 0.05]) if len(vartuple) > 1 and step == ilast: # legends for the additionnal profiles axes = plt.gca() rangex = axes.get_xlim() rangey = axes.get_ylim() xlgd1 = rangex[1] - 0.12 * (rangex[1] - rangex[0]) xlgd2 = rangex[1] - 0.05 * (rangex[1] - rangex[0]) for i in range(profiles.shape[1]): ylgd = rangey[1] - 0.05 * (i + 1) * (rangey[1] - rangey[0]) plt.plot([xlgd1, xlgd2], [ylgd, ylgd], c='black', linestyle=linestyles[i], linewidth=lwdth, dash_capstyle='round',) plt.text(xlgd1 - 0.02 * (rangex[1] - rangex[0]), ylgd, quant[i + 1], ha='right') plt.xlabel(quant[0], fontsize=ftsz) plt.ylabel('z', fontsize=ftsz) plt.xticks(fontsize=ftsz) plt.yticks(fontsize=ftsz) if quant[0] == 'Grid': plt.savefig("Grid" + timename + ".pdf", format='PDF') else: # legend lgd = plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, borderaxespad=0., mode="expand", ncol=3, fontsize=ftsz, columnspacing=1.0, labelspacing=0.0, handletextpad=0.1, handlelength=1.5, fancybox=True, shadow=False) plt.savefig(quant[0].replace(' ', '_') + timename + ".pdf", format='PDF', bbox_extra_artists=(lgd, ), bbox_inches='tight') plt.close(fig) if args.plot_difference: # plot time series of difference profiles if quant[0] == 'Concentration': iminc = concdif.index(min(concdif)) axax[0].semilogy(tsteps[0:ilast:istep, 2], concdif / concdif[0]) axax[0].semilogy(tsteps[iminc * istep, 2], concdif[iminc] / concdif[0], 'o', label=r'$t=%.2e$' % (tsteps[iminc, 2])) axax[0].set_ylabel('Composition diff.') plt.legend(loc='upper right') return tsteps[iminc * istep, 2], concdif[iminc] / concdif[0],\ iminc, timename if quant[0] == 'Temperature': axax[1].semilogy(tsteps[istart:ilast:istep, 2], tempdif / tempdif[0]) imint = tempdif.index(min(tempdif)) axax[1].semilogy(tsteps[imint * istep, 2], tempdif[imint] / tempdif[0], 'o', label=r'$t=%.2e$' % (tsteps[imint, 2])) axax[1].set_ylabel('Temperature diff.') plt.legend(loc='lower right') # maximum velocity as function of time axax[2].semilogy(tsteps[istart:ilast:istep, 2], wmax) axax[2].set_ylabel('Max. rms vert. velocity') axax[2].set_xlabel('Time') wma = max(wmax) iwm = wmax.index(wma) sigma = math.log(wmax[iwm - 3] / wmax[0]) / tsteps[iwm - 3, 2] expw = [wmax[0] * math.exp(sigma * t) for t in tsteps[0:iwm + 2:istep, 2]] axax[2].semilogy(tsteps[0:iwm + 2:istep, 2], expw, linestyle='--', label=r'$sigma=%.2e$' % sigma) plt.legend(loc='upper right') return tsteps[imint * istep, 2], tempdif[imint] / tempdif[0], iwm,\ wma, imint, sigma, timename return None def rprof_cmd(args): """Plot radial profiles""" if not (args.plot_conctheo and args.plot_temperature and args.plot_concentration): args.plot_difference = False ctheoarg = None, None if args.plot_difference: # plot time series of difference profiles # initialize the plot here figd, axax = args.plt.subplots(3, sharex=True) ra0 = args.par_nml['refstate']['ra0'] ctheoarg = axax, None # parameters for the theoretical composition profiles spherical = args.par_nml['geometry']['shape'].lower() == 'spherical' if spherical: rcmb = args.par_nml['geometry']['r_cmb'] else: rcmb = 0. rmin = rcmb rmax = rcmb + 1. rbounds = rmin, rmax, rcmb if args.plot_conctheo: xieut = args.par_nml['tracersin']['fe_eut'] k_fe = args.par_nml['tracersin']['k_fe'] xi0l = args.par_nml['tracersin']['fe_cont'] xi0s = k_fe * xi0l xired = xi0l / xieut rsup = (rmax**3 - xired**(1 / (1 - k_fe)) * (rmax**3 - rmin**3))**(1 / 3) print('rmin, rmax, rsup=', rmin, rmax, rsup) def initprof(rpos): """Theoretical initial profile.""" if rpos < rsup: return xi0s * ((rmax**3 - rmin**3) / (rmax**3 - rpos**3))**(1 - k_fe) else: return xieut ctheoarg = ctheoarg[0], initprof rprof_data = RprofData(args) data, tsteps, nzi = rprof_data.data, rprof_data.tsteps, rprof_data.nzi for var in 'tvnc': # temperature, velo, visco, concentration meta = constants.RPROF_VAR_LIST[var] if not misc.get_arg(args, meta.arg): continue labels = [meta.name] cols = [meta.prof_idx] if misc.get_arg(args, meta.min_max): labels.extend(['Mean', 'Minimum', 'Maximum']) cols.extend([meta.prof_idx + 1, meta.prof_idx + 2]) out = plotprofiles(labels, cols, data, tsteps, nzi, rbounds, args, ctheoarg) if var == 't' and args.plot_difference: _, _, _, _, imint, sigma, timename = out if var == 'c' and args.plot_difference: _, _, iminc, timename = out if args.plot_difference: args.plt.ticklabel_format(style='sci', axis='x') args.plt.savefig('Difference_to_overturned{}.pdf'.format(timename), format='PDF') args.plt.close(figd) with open('statmin.dat', 'w') as fich: fmt = '{:12}' * 6 + '\n' fich.write(fmt.format('rcmb', 'k_fe', 'ra', 'tminT', 'sigma', 'tminC')) fmt = '{:12.5e}' * 6 fich.write(fmt.format(rcmb, k_fe, ra0, tsteps[imint, 2], sigma, tsteps[iminc, 2])) # Plot grid spacing if args.plot_grid: plotprofiles(['Grid'], None, data, tsteps, nzi, rbounds, args, ctheoarg) # Plot the profiles of vertical advection: total and contributions from up- # and down-welling currents if args.plot_advection: plotprofiles(['Advection per unit surface', 'Total', 'down-welling', 'Up-welling'], (57, 58, 59), data, tsteps, nzi, rbounds, args, ctheoarg) if spherical: plotprofiles(['Total scaled advection', 'Total', 'down-welling', 'Up-welling'], (57, 58, 59), data, tsteps, nzi, rbounds, args, ctheoarg, integrate=True) if args.plot_energy: plotprofiles(['Energy', 'Total', 'Advection', 'conduction'], (57, 58, 59), data, tsteps, nzi, rbounds, args, ctheoarg, integrate=True) PKmCH 2DDstagpy/parfile.py"""StagYY par file handling""" import f90nml import os.path from .constants import CONFIG_DIR PAR_DFLT_FILE = os.path.join(CONFIG_DIR, 'par') PAR_DEFAULT = { 'switches': { 'verbose': False, 'try_continue': False, 'compressible': False, 'dens_strat': True, 'dens_actual': False, 'melting': False, 'composition': False, 'phase_changes': False, 'geoid_kernel_mode': False, 'tracers': False, 'cont_rafts': False, 'cont_tracers': False, 'clb_experiment': False, 'chem_bm1': False, 'pvk_bm': False, 'Io_tidal_heating': False, 'multig_solve': True, 'plates_analyse': False, 'dimensional_units': False, 'total_pressure': False, 'dens_read_Fabio': False, 'visc_read_Fabio': False, 'quiet': False, }, 'geometry': { 'shape': 'spherical', 'nxtot': 1, 'nytot': 512, 'nztot': 64, 'aspect_ratio': [1.0, 10.0, 1.0], 'theta_position': 'default', 'zspacing_mode': 'constant', 'dresl_topbot': 1.0, 'D_dimensional': 2890.e3, 'D_dim_top': 0.0e3, 'anti_squeeze': True, 'anti_squeeze_r_fraction': 0.5, 'dresl_air': 1.0, 'dresl_surf': 1.0, 'dresl_cmb': 1.0, 'dresl_ppv': 1.0, 'dresl_660': 1.0, 'dresl_400': 1.0, 'wresl_air': 0.1, 'wresl_surf': 0.1, 'wresl_cmb': 0.1, 'wresl_ppv': 0.1, 'wresl_660': 0.1, 'wresl_400': 0.1, 'r_cmb': 1.19, 'yy_MinOverlap': True, 'kmaxb': 0.25, 'kmaxt': 0.25, 'npbl': 10, }, 'refstate': { 'Di0': 1.18, 'ra0': 1.e5, 'nradiogenic': 1, 'Rh': [0.0], 'Rh_tstart': 0., 'drho_thermal': 0.125, 'perplex_properties': False, 'perplex_hz_file': 'harzburgite.dat', 'perplex_bs_file': 'basalt.dat', 'perplex_pr_file': 'primordial.dat', 'perplex_py_file': 'pyrolite.dat', 'perplex_dT': 100., 'Tref_surf': 0.64, 'rhoref_surf': [3240., 3080., 2900., 2900., 3080., 3080., 7000.], 'gr_s': [[1.3]], 'drho_surf': 1.0, 'drho_cmb': 4337., 'q_gamma': [[1.0]], 'delT0_expan': [[6.0]], 'ex_tkappa': [[3.0]], 'ak_expan': [1.4], 'c_dependent_heating': False, 'dHeat_dprim': 1.0, 'dHeat_dbasalt': 1.0, 'Rh_halflife': ['BIG'], 'deltaT_dimensional': 2500.0, 'g_dimensional': 9.81, 'gravity_3D': False, 'core_superheat': 0.0, 'tcond_dimensional': 3.0, 'core_tcond': 46.0, 'dens_dimensional': 3300., 'expan_dimensional': 3.e-5, 'Cp_dimensional': 1200., 'core_cooling': False, 'core_model': 'simple', 'core_Kppm': 0.0, 'core_Tmelt0': 5600, 'core_heatcap': 6.38e12, 'phase_assemblage_varies_with_C': False, 'air_density': 1.e-3, 'topregion_adiabatic': False, 'topregion_z': 1.0, 'topregion_thalf': 'BIG', 'Birch_Murnaghan': False, 'K0': [[210.e9]], 'Kp': [[3.9]], 'heterogeneous_heating': False, 'Hs_zbot': 0., 'Hs_ztop': 1., 'Hs_bot': 0., 'Hs_top': 0., 'Hs_mode': 1, }, 'boundaries': { 'topT_mode': 'isothermal', 'topT_val': 1.0, 'botT_mode': 'iso', 'botT_val': 0.0, 'outT_val': -1.0, 'topV_mode': 'free-slip', 'botV_mode': 'free-slip', 'x_patch': 0.5, 'y_patch': 0.5, 'r_patch': 0.1, 'x_spread': 0.0, 'v_spread': 0.0, 'nx_plates': 2, 'ny_plates': 2, 'x_bc': 'ww', 'y_bc': 'ww', 'v_mantle': 0.0, 'air_layer': False, 'air_thickness': 0.1, 'air_threshold': 0.8, 'virtual_bndry_distance': -1., 'zero_surface_mean_flow': True, 'Tsurf_eqm': 300., 'topT_locked_subsolar': 300., 'topT_locked_farside': 200., 'vbc_file_stem': 'vbc', 'plateid_file_stem': 'pid', 'vbc_file_interval_Myr': 1.0, 'vbc_file_start_Myr': -70., 'platerecon_time_origin': 0.0, 'platerecon_time_init': 0.0, 'platerecon_time_end': 0.0, 'BottomPhaseChange': False, 'BotPphase': 1., 'TopPhaseChange': False, 'TopPphase': 1., }, 't_init': { 'imode_t': 'cells', 'amp_t': 0.00, 'blthick': 0.1, 't0_init': 0.5, 'z0_t': 0.5, 'w0_t': 0.5, 'dip': 0.0, 'length': 0.0, 'ocean_age': 60.0, 'cont_thick': 100.e3, }, 'timein': { 'nsteps': 1, 'nwrite': 1, 'alpha_adv': 0.5, 'alpha_diff': 1.0, 'iord': 2, 'advection_scheme': 'MPDATA', 'stoptime': 'hhmm', 'dt_write': 'BIG', 'finish_at_time': 'BIG', 'TVD_limiter': 2., 'diffusion_implicit': False, 'diffusion_beta': 1.0, 'diffusion_errmax': 1.e-4, 'diffusion_maxits': 100, 'diffusion_alpha': 0.6, 'diffusion_nrelax': 3, 'density_jump_stabilization': False, 'density_jump_theta': 1.0, 'density_jump_simplify': True, 'PBSTimeRequested': 168.0, }, 'viscosity': { 'ietalaw': 0, 'E_eta': [[[0.0]]], 'V_eta': [[[0.0]]], 'Pdecay_V': [[['BIG']]], 'Ts_eta': -1.0, 'n_eta': 3.5, 'eta0': 1.0, 'T0_eta': 0.5, 'd0_eta': 0.5, 'stress0_eta': [[1.0]], 'deta_p': [[1.0]], 'deta_p_correct_jumps': False, 'diffusion_creep': True, 'dislocation_creep': False, 'eta_max': 1.e5, 'eta_min': 1.e-5, 'eta_minmax_depthmax': 'BIG', 'deta_basalt': 1.0, 'deta_ccrust': 1.0, 'deta_primordial': 1.0, 'deta_metal': 1.0, 'deta_water': 1.0, 'etalaw_basalt': 1, 'etalaw_ccrust': 1, 'etalaw_metal': 1, 'etalaw_melt': 1, 'etalaw_primordial': 1, 'etalaw_water': 1, 'c0eta_basalt': 0.0, 'c0eta_ccrust': 0.0, 'c0eta_primordial': 0.0, 'c0eta_metal': 0.0, 'c0eta_water': 0.0, 'weak_line': False, 'weak_dx': 0.1, 'weak_dz': 0.05, 'weak_eta': 1.0, 'strong_cont_rafts': False, 'strong_eta': 1.0, 'strong_dz': 0.05, 'vertonly_eta': False, 'horzonly_eta': False, 'plasticity': False, 'c_dependent_yieldpars': False, 'stressY_eta': 'BIG', 'dstressY_eta_dp': 0.0, 'cohesion_eta': 'BIG', 'frictionCoeff_eta': 'BIG', 'stressY_basalt': 'BIG', 'dstressY_basalt': 0.0, 'cohesion_basalt': 'BIG', 'frictionCoeff_basalt': 'BIG', 'stressY_ccrust': 'BIG', 'dstressY_ccrust': 0.0, 'cohesion_ccrust': 'BIG', 'frictionCoeff_ccrust': 'BIG', 'stressY_metal': 'BIG', 'dstressY_metal': 0.0, 'cohesion_metal': 'BIG', 'frictionCoeff_metal': 'BIG', 'stressY_primordial': 'BIG', 'dstressY_primordial': 0.0, 'cohesion_primordial': 'BIG', 'frictionCoeff_primordial': 'BIG', 'd_etalaw': 1, 'Rheal_law': 'constant', 'E_Rh': 0.0, 'F_eta': 0.0, 'eta_melt': False, 'deta_melt': 1.0, 'mean_T_eta': False, 'T0_Tmean': 0.5, 'E_eta_Tmean': 0.0, 'eta_conv': 0.01, 'alpha_eta': 1.0, 'maxits_eta': 20, 'bmoore_eta': False, 'aniso_ep': False, 'visc_interp_law': 'g', 'phase_dep_Eeta': False, 'eta_air': 1., 'PressureDepRheology': False, 'PressureDepVact': False, 'RaFromRheolPref': False, }, 'grainsize': { 'graingrowth_k': 1., 'graingrowth_p': 4., 'graingrowth_E': 300.e3, 'grainreduction_C': 1., 'grainsize0_eta': 1., 'grainsize_n_eta': 3.0, 'GS_ForcedPiezometric': True, 'GS_TwoPhases': False, 'GS_PiezometerModel': 'homogeneous', 'GG_k0': [2.49813e7, 2.49813e7], 'GG_p': [2.0, 2.0], 'GG_E': [2.e5, 2.e5], 'GS_InitSize': [1e2, 1e2], 'GS_gamma': [1e6, 1e6], 'curv_k0': 74943.9, 'curv_E': 2.e5, 'curv_q': 1.5, 'curv_gamma': 1.e6, 'curv_Init': 1e2, 'GSD_sigma': 0.8, 'GS_fG0': 1.e-2, 'GS_fI0': 0.0, 'transition_depth_GS': [410.e3, 500.e3, 660.e3, 2740.e3], 'GSAfterPhaseTrans': 5.e0, 'curvDownThrough660': 5.e0, }, 'iteration': { 'maxits': 50, 'minits': 1, 'errmax': 1.e-3, 'alpha_mom': 0.7, 'alpha_cont': 0.9, 'redblack': False, 'extra_lid_relax': False, 'nlid_relax': 1, 'tlid': 0.2, 'extra_hiresid_relax': False, 'nhi_relax': 1, 'normalize_res_by_eta': False, 'relax_kernel': 'classic', }, 'multi': { 'maxcoarseits': 200, 'nrelax1': 3, 'nrelax2': 3, 'nhmin': 4, 'nzmin': 4, 'nzpnmin': 1, 'alpha_multi': 1.0, 'matrix_dep_pinterp': True, 'matrix_dep_vinterp': False, 'ncpus_per_node': 8, 'nppcpu_opt_samenode': 64, 'nppcpu_opt_xnode': 1024, 'more_coarse_its': True, 'f_cycle': False, 'direct_coarse_soln': False, 'idle_duplicates': True, }, 'ioin': { 'input_file': 'null', 'input_frame': 1, 'output_file_stem': 'test', 't_write': True, 'vp_write': False, 'eta_write': False, 'c_write': False, 'bs_write': False, 'hz_write': False, 'air_write': False, 'cc_write': False, 'primordial_write': False, 'f_write': False, 't_spectrum_write': False, 'vm_write': False, 'stress_write': False, 'g_write': False, 'tra_write': False, 'metal_write': False, 'vd_write': False, 'edot_write': False, 'age_write': False, 'nmelt_write': False, 'rho_write': False, 'grainsize_write': False, 'hpe_write': False, 'water_write': False, 'K_Ar_write': False, 'oxides_write': False, 'ph_write': False, 'residue_write': False, 'nrho_write': False, 'w_write': False, 'delete_old_trafiles': True, 'delete_old_nonhdf_files': False, 'torpol_write': False, 'svel_write': False, 'pvel_write': False, 'bvel_write': False, 'topo_write': False, 'lyap_write': False, 'divvor_write': False, 'crust_write': False, 'heatflux_write': False, 'heatflux_spectrum_write': False, 'lmax': 16, 'stress_axis_write': False, 'trapercell_write': False, 'write_32bit_precision': False, 'MORBprof_write': False, 'save_file_framestep': 10, 'overwrite_old_files': False, 'restart_from_meshed_fields_only': False, 'hdf5_output_folder': '+hdf5', 'hdf5_input_folder': '+hdf5', 'MaxFileSizeAllowed': 100000000, 'MaxFileSizeAllowedTra': 100000000, }, 'compin': { 'imode_comp': 'layered', 'zcomp': 0.5, 'amp_comp': 0.0, 'ra_comp': [0.0], 'lenardic_filter': False, 'outC_val': -1.0, 'dzcomp': 0.1, }, 'melt': { 'imode_melt': 'whatever', 'ra_melt': 0.0, 'mr_melt': 1.e6, 'Tsol0': 0.5, 'dTsol_dz': 1.0, 'latent_heat': 1.0, 'amp_melt': 0.02, 'eruption': False, 'erupt_fraction': 1.0, 'intrusion': False, 'intrude_fraction': 1.0, 'tzintrusion': False, 'tzintrude_fraction': 1.0, 'outF_val': -1.0, 'erupt_mode': 'all', 'solidus_function': 'simple', 'd_intrude': 0.0, 'intrude_position': 'surface', 'harzburgite_melts': False, 'metal_melts': True, 'suppress_asth': False, 'c_dependent_solidus': False, 'd_erupt_Earth': 1.0, 'erupt_threshold': 0.0, 'deltaTsol_depletion_dimensional': 60., 'deltaTsol_water_dimensional': 0., 'ddens_sol_liq_dimensional': 500., 'eta_liquid_dimensional': 10., 'k0_dimensional': 3.e-10, 'p_dependent_solidus': False, 'Cf_dependent_Tmelting': False, 'melt_phase_systems': False, 'melt_segregation': False, 'erupt_T_lith': 1400., 'magma_effective_tcond': False, 'magma_tcond_factor': 1.e5, 'magma_topTBL_const': 1.e-7, 'fractional_xtallisation_n_melting': False, 'trackTTGsource': False, }, 'phase': { 'nphase_systems': 1, 'Ol_frac_ref': 0.6, 'basalt_frac_ref': 0.2, 'effective_Cp_alpha': True, 'scale_PCdepths': False, 'g0_dimensional_PCdepths': 9.81, 'D0_dimensional_PCdepths': 2890.e3, 'sys_ol': 1, 'sys_pxgt': 2, 'sys_melt_ol': 3, 'sys_melt_pxgt': 4, 'sys_prim': 5, 'sys_ccrust': 6, 'sys_metal': 7, }, 'continents': { 'ncont': 1, 'cont_radius': [0.1], 'cont_startpos': [[0.5]], 'cont_offsetplot': True, 'rafts_move': False, 'cont_collideaction': 'ignore', 'cont_Biot': 1.e6, 'cont_thickness': 100.e3, }, 'tracersin': { 'ntracers': 50, 'ntracers_per_cell': -1, 'imode_tra': 'even+ran', 'B_basalt': 0.0, 'B_ccrust': 0.0, 'B_metal': 0.0, 'd_crust': 0.1, 'w_crust': 0.0, 'Ctr_truncate': False, 'k_fe': 0.85, 'Fe_eut': 0.8, 'Fe_cont': 0.1, 'Ctr_lininterp': True, 'Ttr_lininterp': True, 'trastrain_diagnostics': False, 'tracers_everywhere': True, 'tracers_strain': False, 'tracers_timemelt': False, 'tracers_recstartpos': False, 'tracers_advord': 2, 'tracers_interp2ord': True, 'tracers_weakcrust': False, 'weakcrust_maxdepth': 0.3, 'tracers_temperature': False, 'oxides_compositions': False, 'basalt_harzburgite': False, 'continental_crust': False, 'tracers_weakfault': False, 'traceelement_hpe': False, 'traceelement_water': False, 'traceelement_K_Ar': False, 'track_GrainSize': False, 'metal': False, 'ddens_SiO2': 0.0, 'ddens_MgO': 0.0, 'ddens_FeO': 0.0, 'ddens_XO': 0.0, 'dTmelting_dFeO': 0.0, 'primordial_layer': True, 'd_primordial': 0.1, 'B_primordial': 0.0, 'imode_primordial': 'layer', 'imode_metal': 'top_layer', 'd_metal': 0.0, 'tracers_horizontaladvection': False, 'tracers_noadvection': False, 'Dpartition_FeO': 1.0, 'Dpartition_hpe': 1.0, 'Dpartition_water': 1.0, 'Dpartition_K': 1.0, 'Dpartition_Ar': 1.0, 'outgas_fraction_water': 1.0, 'outgas_fraction_Ar': 1.0, 'ingas_water': False, 'ingas_depth': 10.e3, 'IntrudedBasaltHydrationFraction': 0.0, 'overturn': False, 'mode_init_overturn': 'linear', }, 'plot': { 'plot_file_stem': 'image', 'npix': 200000, 'dots': True, 'auto_t_scale': True, 'z_xyslice': 0.5, 't_superadiabatic': False, 't_plot': True, 'v_plot': False, 'p_plot': False, 'eta_plot': False, 'rho_plot': False, 'c_plot': False, 'f_plot': False, 'pd_plot': False, 'dT_plot': False, 'vm_plot': False, 'stress_plot': False, 'edot_plot': False, 'g_plot': False, 'tra_plot': False, 'vd_plot': False, 'air_plot': False, 'oxides_plot': False, 'hpe_plot': False, 'water_plot': False, 'K_Ar_plot': False, 'primordial_plot': False, 'rhs_plot': False, 'bs_plot': False, 'hz_plot': False, 'cc_plot': False, 'tcond_plot': False, 'Cp_plot': False, 'expan_plot': False, 'grainsize_plot': False, 'age_plot': False, 'nmelt_plot': False, 'strain_plot': False, 'ph_plot': False, 'lyap_plot': False, 'residue_plot': False, 'trapercell_plot': False, 'imgtype': 'png', }, } def _write_default(): """create default par file""" if not os.path.isfile(PAR_DFLT_FILE): f90nml.write(PAR_DEFAULT, PAR_DFLT_FILE) def _read_default(): """read default par file""" _write_default() return f90nml.read(PAR_DFLT_FILE) def readpar(args): """read StagYY par file""" par_file = os.path.join(args.path, 'par') par_dflt = _read_default() if os.path.isfile(par_file): par_nml = f90nml.read(par_file) for section in par_nml: par_dflt[section].update(par_nml[section]) else: if not (args.create or args.update or args.edit): print('no par file found, check path') par_nml = par_dflt return par_nml PK+lHH%stagpy/constants.py"""defines some constants""" import os.path from collections import OrderedDict, namedtuple CONFIG_DIR = os.path.expanduser('~/.config/stagpy') Varf = namedtuple('Varf', ['par', 'name', 'arg']) FIELD_VAR_LIST = OrderedDict(( ('t', Varf('t', 'Temperature', 'plot_temperature')), ('c', Varf('c', 'Composition', 'plot_composition')), ('n', Varf('eta', 'Viscosity', 'plot_viscosity')), ('d', Varf('rho', 'Density', 'plot_density')), ('h', Varf('wtr', 'Water', 'plot_water')), ('a', Varf('age', 'Age', 'plot_age')), ('u', Varf('vp', 'x Velocity', 'plot_xvelo')), ('v', Varf('vp', 'y Velocity', 'plot_yvelo')), ('w', Varf('vp', 'z Velocity', 'plot_zvelo')), ('p', Varf('vp', 'Pressure', 'plot_pressure')), ('s', Varf('vp', 'Stream function', 'plot_stream')), )) Varr = namedtuple('Varr', ['name', 'arg', 'min_max', 'prof_idx']) RPROF_VAR_LIST = OrderedDict(( ('t', Varr('Temperature', 'plot_temperature', 'plot_minmaxtemp', 1)), ('v', Varr('Vertical velocity', 'plot_velocity', 'plot_minmaxvelo', 7)), ('n', Varr('Viscosity', 'plot_viscosity', 'plot_minmaxvisco', 13)), ('c', Varr('Concentration', 'plot_concentration', 'plot_minmaxcon', 36)), ('g', Varr('Grid', 'plot_grid', None, None)), ('a', Varr('Advection', 'plot_advection', None, None)), ('e', Varr('Energy', 'plot_energy', None, None)), ('h', Varr('Concentration Theo', 'plot_conctheo', None, None)), ('i', Varr('Init overturn', 'plot_overturn_init', None, None)), ('d', Varr('Difference', 'plot_difference', None, None)), )) PKXILH~ʧ stagpy/misc.py"""miscellaneous definitions""" import importlib from itertools import zip_longest from math import ceil import os.path import sys INT_FMT = '{:05d}' def stop(*msgs): """print error message and exit""" print('ERROR:', *msgs, file=sys.stderr) sys.exit() def _file_name(args, fname): """return full name of StagYY out file""" return os.path.join(args.path, args.name + '_' + fname) def stag_file(args, fname, timestep=None, suffix=''): """return name of StagYY out file if exists specify a time step if needed """ if timestep is not None: fname = fname + INT_FMT.format(timestep) fname = _file_name(args, fname + suffix) if not os.path.isfile(fname): stop('requested file {} not found'.format(fname)) return fname def out_name(args, par_type): """return out file name format for any time step""" return args.outname + '_' + par_type + INT_FMT def set_arg(args, arg, val): """set a cmd line with arg string name""" vars(args)[arg] = val def get_arg(args, arg): """set a cmd line with arg string name""" return vars(args)[arg] def lastfile(args, begstep): """look for the last binary file research based on temperature files """ fmt = _file_name(args, 't' + INT_FMT) endstep = 100000 while begstep + 1 < endstep: guess = int(ceil((endstep + begstep) / 2)) if os.path.isfile(fmt.format(guess)): begstep = guess else: endstep = guess return begstep def parse_line(line, convert=None): """convert columns of a text line line values have to be space separated, values are converted to float by default. convert argument is a list of functions used to convert the first values. """ if convert is None: convert = [] line = line.split() for val, func in zip_longest(line, convert[:len(line)], fillvalue=float): yield func(val) def parse_timesteps(args): """parse timestep argument""" tstp = args.timestep.split(':') if not tstp[0]: tstp[0] = '0' if len(tstp) == 1: tstp.extend(tstp) if not tstp[1]: tstp[1] = lastfile(args, int(tstp[0])) tstp[1] = int(tstp[1]) + 1 if len(tstp) != 3: tstp = tstp[0:2] + [1] if not tstp[2]: tstp[2] = 1 args.timestep = list(map(int, tstp)) def plot_backend(args): """import matplotlib and seaborn""" mpl = importlib.import_module('matplotlib') if args.matplotback: mpl.use(args.matplotback) plt = importlib.import_module('matplotlib.pyplot') if args.useseaborn: sns = importlib.import_module('seaborn') args.sns = sns if args.xkcd: plt.xkcd() args.mpl = mpl args.plt = plt PKmCHMstagpy/stagpy.py""" Read and plot stagyy binary data Author: Martina Ulvrova Date: 2014/12/02 """ from . import config def main(): """stagpy entry point""" args = config.parse_args() args.func(args) PKQLH&stagpy-0.1.0.dist-info/DESCRIPTION.rst.. image:: https://landscape.io/github/mulvrova/StagPy/master/landscape.svg?style=flat-square :target: https://landscape.io/github/mulvrova/StagPy/master :alt: Code Health StagPy ====== StagPy is a Python 3 command line tool to read and process StagYY output files to produce high-quality figures. The aim is to have different cases in one file (Cartesian, Spherical Annulus, etc). The code to read the binary output files has been adapted from a matlab version initially developed by Boris Kaus. Installation ============ *if you want to use (and modify) the development version, see the "For developers" section at the end of this page* You will need Python 3.3 or higher to use StagPy. If you don't have ``pip`` for Python3 on your system, download the official script and run it with ``python3``. StagPy is available via ``pip``. You can install it with the following command:: python3 -m pip install --user stagpy Make sure that the directory where ``pip`` install package entry-points (usually ``~/.local/bin``) is in your ``PATH`` environment variable. You can run ``python3 -m pip show stagpy`` to obtain some hint about this location (this command will show you were the compiled sources are installed, e.g. ``~/.local/lib/python3.5/site-packages``, from which you can deduce the entry-point location, e.g. ``~/.local/bin``). Once this is done, you can enable command-line auto-completion if you use either bash or zsh. Add this to your ``~/.bashrc`` file:: eval $(register-python-argcomplete stagpy) Or this to your ``~/.zshrc`` file:: autoload bashcompinit bashcompinit eval $(register-python-argcomplete stagpy) Finally, run the following once to create your config file (at ``~/.config/stagpy/``):: stagpy config --create Enjoy! Available commands ================== The available subcommands are the following: * ``field``: computes and/or plots scalar fields such as temperature or stream function; * ``rprof``: computes and/or plots radial profiles; * ``time``: computes and/or plots time series; * ``plates``: plate analysis; * ``var``: displays a list of available variables; * ``config``: configuration handling. You can run ``stagpy --help`` (or ``stagpy -h``) to display a help message describing those subcommands. You can also run ``stagpy --help`` to have some help on the available options for one particular sub command. StagPy looks for a StagYY ``par`` file in the current directory. It then reads the value of the ``output_file_stem`` option to determine the location and name of the StagYY output files (set to ``test`` if no ``par`` file can be found). You can change the directory in which StagYY looks for a ``par`` file by two different ways: * you can change the default behavior in a global way by editing the config file (``stagpy config --edit``) and change the ``core.path`` variable; * or you can change the path only for the current run with the ``-p`` option. The time step option ``-s`` allows you to specify a range of time steps in a way which mimic the slicing syntax: ``begin:end:gap`` (both ends included). If the first step is not specified, it is set to ``0``. If the final step is not specified, all available time steps are processed. Here are some examples: * ``-s 100:350`` will process every time steps between 100 and 350; * ``-s 201:206:2`` will process time steps 201, 203 and 205; * ``-s 201:205:2`` same as previous; * ``-s 5682:`` will process every time steps from the 5682nd to the last one; * ``-s :453`` will process every time steps from the 0th to the 453rd one; * ``-s ::2`` will process every even time steps. By default, the temperature, pressure and stream function fields are plotted. You can change this with the ``-o`` option (e.g. ``./main.py field -o ps`` to plot only the pressure and stream function fields). For developers ============== A ``Makefile`` in the git repository allows you to install StagPy in a virtual environment. StagPy uses the following non-standard modules: numpy, scipy, f90nml, matplotlib, and seaborn (the latter is optional and can be turned off with the ``core.useseaborn`` option). These dependencies will be checked and needed installation performed in a virtual environment. If you use Python3.2 or encouter problems with the installation, see the troubleshooting section at the end of this README. However, installation of ``numpy`` and ``scipy`` involve heavy building operations, it might be better that you (or your system administrator) install it with a package manager such as ``homebrew`` on Mac OS or your favorite Linux package manager. The installation process is then fairly simple:: git clone https://github.com/mulvrova/StagPy.git cd StagPy make A soft link named ``stagpy`` is created in your ``~/bin`` directory, allowing you to launch StagPy directly by running ``stagpy`` in a terminal (provided that ``~/bin`` is in your ``PATH`` environment variable). Two files ``.comp.zsh`` and ``.comp.sh`` are created. Source them respectively in ``~/.zshrc`` and ``~/.bashrc`` to enjoy command line completion with zsh and bash. Run ``make info`` to obtain the right sourcing commands. To check that everything work fine, go to the ``data`` directory of the repository and run:: stagpy field Three PDF files with a plot of the temperature, pressure and stream function fields should appear. Troubleshooting =============== * Matplotlib related error in MacOS This might be due to the matplotlib backend that is not correctly set. See this Stack Overflow question: * Installation fails with ``ImportError: No module named 'encodings'`` This seems to be due to a bug in the venv module with some Python installation setups. If installing Python properly with your package manager doesn't solve the issue, you can try installing StagPy without any virtual environment by using ``make novirtualenv``. PKQLHFKI//'stagpy-0.1.0.dist-info/entry_points.txt[console_scripts] stagpy = stagpy.stagpy:main PKQLH -$stagpy-0.1.0.dist-info/metadata.json{"classifiers": ["Development Status :: 3 - Alpha", "Intended Audience :: Science/Research", "License :: OSI Approved :: GNU General Public License v2 (GPLv2)", "Programming Language :: Python :: 3 :: Only", "Programming Language :: Python :: 3.3", "Programming Language :: Python :: 3.4", "Programming Language :: Python :: 3.5"], "extensions": {"python.commands": {"wrap_console": {"stagpy": "stagpy.stagpy:main"}}, "python.details": {"contacts": [{"name": "Martina Ulvrova, Adrien Morison, St\u00e9phane Labrosse", "role": "author"}], "document_names": {"description": "DESCRIPTION.rst"}, "project_urls": {"Home": "https://github.com/mulvrova/StagPy"}}, "python.exports": {"console_scripts": {"stagpy": "stagpy.stagpy:main"}}}, "extras": [], "generator": "bdist_wheel (0.26.0)", "license": "GPLv2", "metadata_version": "2.0", "name": "stagpy", "run_requires": [{"requires": ["argcomplete", "f90nml", "matplotlib", "numpy", "scipy", "seaborn"]}], "summary": "Tool for StagYY output files processing", "version": "0.1.0"}PKQLH%$stagpy-0.1.0.dist-info/top_level.txtstagpy PKQLH}\\stagpy-0.1.0.dist-info/WHEELWheel-Version: 1.0 Generator: bdist_wheel (0.26.0) Root-Is-Purelib: true Tag: py3-none-any PKQLHr{stagpy-0.1.0.dist-info/METADATAMetadata-Version: 2.0 Name: stagpy Version: 0.1.0 Summary: Tool for StagYY output files processing Home-page: https://github.com/mulvrova/StagPy Author: Martina Ulvrova, Adrien Morison, Stéphane Labrosse Author-email: UNKNOWN License: GPLv2 Platform: UNKNOWN Classifier: Development Status :: 3 - Alpha Classifier: Intended Audience :: Science/Research Classifier: License :: OSI Approved :: GNU General Public License v2 (GPLv2) Classifier: Programming Language :: Python :: 3 :: Only Classifier: Programming Language :: Python :: 3.3 Classifier: Programming Language :: Python :: 3.4 Classifier: Programming Language :: Python :: 3.5 Requires-Dist: argcomplete Requires-Dist: f90nml Requires-Dist: matplotlib Requires-Dist: numpy Requires-Dist: scipy Requires-Dist: seaborn .. image:: https://landscape.io/github/mulvrova/StagPy/master/landscape.svg?style=flat-square :target: https://landscape.io/github/mulvrova/StagPy/master :alt: Code Health StagPy ====== StagPy is a Python 3 command line tool to read and process StagYY output files to produce high-quality figures. The aim is to have different cases in one file (Cartesian, Spherical Annulus, etc). The code to read the binary output files has been adapted from a matlab version initially developed by Boris Kaus. Installation ============ *if you want to use (and modify) the development version, see the "For developers" section at the end of this page* You will need Python 3.3 or higher to use StagPy. If you don't have ``pip`` for Python3 on your system, download the official script and run it with ``python3``. StagPy is available via ``pip``. You can install it with the following command:: python3 -m pip install --user stagpy Make sure that the directory where ``pip`` install package entry-points (usually ``~/.local/bin``) is in your ``PATH`` environment variable. You can run ``python3 -m pip show stagpy`` to obtain some hint about this location (this command will show you were the compiled sources are installed, e.g. ``~/.local/lib/python3.5/site-packages``, from which you can deduce the entry-point location, e.g. ``~/.local/bin``). Once this is done, you can enable command-line auto-completion if you use either bash or zsh. Add this to your ``~/.bashrc`` file:: eval $(register-python-argcomplete stagpy) Or this to your ``~/.zshrc`` file:: autoload bashcompinit bashcompinit eval $(register-python-argcomplete stagpy) Finally, run the following once to create your config file (at ``~/.config/stagpy/``):: stagpy config --create Enjoy! Available commands ================== The available subcommands are the following: * ``field``: computes and/or plots scalar fields such as temperature or stream function; * ``rprof``: computes and/or plots radial profiles; * ``time``: computes and/or plots time series; * ``plates``: plate analysis; * ``var``: displays a list of available variables; * ``config``: configuration handling. You can run ``stagpy --help`` (or ``stagpy -h``) to display a help message describing those subcommands. You can also run ``stagpy --help`` to have some help on the available options for one particular sub command. StagPy looks for a StagYY ``par`` file in the current directory. It then reads the value of the ``output_file_stem`` option to determine the location and name of the StagYY output files (set to ``test`` if no ``par`` file can be found). You can change the directory in which StagYY looks for a ``par`` file by two different ways: * you can change the default behavior in a global way by editing the config file (``stagpy config --edit``) and change the ``core.path`` variable; * or you can change the path only for the current run with the ``-p`` option. The time step option ``-s`` allows you to specify a range of time steps in a way which mimic the slicing syntax: ``begin:end:gap`` (both ends included). If the first step is not specified, it is set to ``0``. If the final step is not specified, all available time steps are processed. Here are some examples: * ``-s 100:350`` will process every time steps between 100 and 350; * ``-s 201:206:2`` will process time steps 201, 203 and 205; * ``-s 201:205:2`` same as previous; * ``-s 5682:`` will process every time steps from the 5682nd to the last one; * ``-s :453`` will process every time steps from the 0th to the 453rd one; * ``-s ::2`` will process every even time steps. By default, the temperature, pressure and stream function fields are plotted. You can change this with the ``-o`` option (e.g. ``./main.py field -o ps`` to plot only the pressure and stream function fields). For developers ============== A ``Makefile`` in the git repository allows you to install StagPy in a virtual environment. StagPy uses the following non-standard modules: numpy, scipy, f90nml, matplotlib, and seaborn (the latter is optional and can be turned off with the ``core.useseaborn`` option). These dependencies will be checked and needed installation performed in a virtual environment. If you use Python3.2 or encouter problems with the installation, see the troubleshooting section at the end of this README. However, installation of ``numpy`` and ``scipy`` involve heavy building operations, it might be better that you (or your system administrator) install it with a package manager such as ``homebrew`` on Mac OS or your favorite Linux package manager. The installation process is then fairly simple:: git clone https://github.com/mulvrova/StagPy.git cd StagPy make A soft link named ``stagpy`` is created in your ``~/bin`` directory, allowing you to launch StagPy directly by running ``stagpy`` in a terminal (provided that ``~/bin`` is in your ``PATH`` environment variable). Two files ``.comp.zsh`` and ``.comp.sh`` are created. Source them respectively in ``~/.zshrc`` and ``~/.bashrc`` to enjoy command line completion with zsh and bash. Run ``make info`` to obtain the right sourcing commands. To check that everything work fine, go to the ``data`` directory of the repository and run:: stagpy field Three PDF files with a plot of the temperature, pressure and stream function fields should appear. Troubleshooting =============== * Matplotlib related error in MacOS This might be due to the matplotlib backend that is not correctly set. See this Stack Overflow question: * Installation fails with ``ImportError: No module named 'encodings'`` This seems to be due to a bug in the venv module with some Python installation setups. If installing Python properly with your package manager doesn't solve the issue, you can try installing StagPy without any virtual environment by using ``make novirtualenv``. PKQLH8!stagpy-0.1.0.dist-info/RECORDstagpy/__init__.py,sha256=47DEQpj8HBSa-_TImW-5JCeuQeRkm5NMpJWZG3hSuFU,0 stagpy/__main__.py,sha256=qGEWJ0bbc1Hmt2OIXR9bSbZTNWpPkrW3TLdli87YuYg,87 stagpy/commands.py,sha256=HkM1UpWP8HgeBxodckq1PU8N5KylWb2z5jOr13m7Xis,1655 stagpy/config.py,sha256=GBJtYtHUEBOi8dToflZsQnY0DPnNIaq6ms4-yvFXMg4,13178 stagpy/constants.py,sha256=5DxanbCYz8g-D65CGL_mU5Hz3EZQYB44sWsrXRelxVo,1567 stagpy/field.py,sha256=NB7PuzZwVNrs1Vp4_eQbL1RMtzQGOlHlourxsDb4pm0,2713 stagpy/misc.py,sha256=u32JzDi2jl5-MU4c1wF2MQ5SJgKAFSLJcPrI2ZLhehM,2777 stagpy/parfile.py,sha256=V44B--wTZDslMfIz5VWBTLnbcvAzgCn_0eqfM9u3Vj4,17566 stagpy/plates.py,sha256=vIbfmH_SzDmoxiE4aBY5TgdjDk4YTwibFGNjJHJ5UIM,18741 stagpy/rprof.py,sha256=DaZmd_F9EKbgMplpkWEyEbPPbdugwpmAr0FBH8mUEoY,15111 stagpy/stagdata.py,sha256=VpXyfClPfzMYD3yZHbyyJOSk6JnPcOWzJhYG9HoL4Qw,10149 stagpy/stagpy.py,sha256=LttYuGYyXtbEPjFCin1kjDaPGzgEkFtRIf0O94BzxGA,199 stagpy/time_series.py,sha256=R7ArRXd64Hq45fQ8YipWWC1ZK7Rhf5-JuKfCBp73Pcc,4602 stagpy-0.1.0.dist-info/DESCRIPTION.rst,sha256=-zP-8IjcO-urXeT0Qc-tBbhK5ItwsZYVbcS47qDd0c0,6085 stagpy-0.1.0.dist-info/METADATA,sha256=urxYyNW3YRvvZG61CYu8BDwd79TX4pKBOh8NDX8tkM0,6863 stagpy-0.1.0.dist-info/RECORD,, stagpy-0.1.0.dist-info/WHEEL,sha256=zX7PHtH_7K-lEzyK75et0UBa3Bj8egCBMXe1M4gc6SU,92 stagpy-0.1.0.dist-info/entry_points.txt,sha256=SdT8j2EP97ArdqAjbwKU2njL16Ihz7eNv6yDej9b4JM,47 stagpy-0.1.0.dist-info/metadata.json,sha256=WtILMAgxcP5WHXNT1MKOv8jtiaWAiF94YS1SJMIZHNM,1022 stagpy-0.1.0.dist-info/top_level.txt,sha256=bqMfpWCVzAAEx1Sw__lU_WeBuoDx1LwVK7x181VISOM,7 PKrCHXI{wwstagpy/commands.pyPKmFHHp= stagpy/time_series.pyPK4HWWstagpy/__main__.pyPKXILHz3z3[stagpy/config.pyPKXILH=G~ Mstagpy/field.pyPKXILHqqhQ''Wstagpy/stagdata.pyPKI3Hstagpy/__init__.pyPKXILHoq5I5Istagpy/plates.pyPKXILHGc;;1stagpy/rprof.pyPKmCH 2DDestagpy/parfile.pyPK+lHH%2Istagpy/constants.pyPKXILH~ʧ Ostagpy/misc.pyPKmCHMZstagpy/stagpy.pyPKQLH&|[stagpy-0.1.0.dist-info/DESCRIPTION.rstPKQLHFKI//'sstagpy-0.1.0.dist-info/entry_points.txtPKQLH -$sstagpy-0.1.0.dist-info/metadata.jsonPKQLH%$9xstagpy-0.1.0.dist-info/top_level.txtPKQLH}\\xstagpy-0.1.0.dist-info/WHEELPKQLHr{ystagpy-0.1.0.dist-info/METADATAPKQLH8!$stagpy-0.1.0.dist-info/RECORDPKbs