# -*- coding: utf-8 -*-
import glob
import re
import os
import copy
import numpy as np
import matplotlib.pyplot as plt
from .npfile import npfile
from .movie import Movie
from .log import MagicSetup
from .libmagic import scanDir
[docs]class TOMovie(Movie):
"""
This class allows to read and display the :ref:`TO_mov.TAG <secTO_movieFile>`
generated when :ref:`l_TOmovie=.true. <varl_TOmovie>` is True.
>>> # This will allow you to pick up one TO_mov files among the existing ones
>>> t = TOMovie()
>>> # Read TO_mov.N0m2, time-averaged it and display it with 65 contour levels
>>> t = TOMovie(file='TO_mov.N0m2', avg=True, levels=65, cm='seismic')
"""
[docs] def __init__(self, file=None, iplot=True, cm='seismic', datadir='.',
cut=0.5, levels=33, avg=True, precision=np.float32,
lastvar=None, nvar='all'):
"""
:param file: the filename of the TO_mov file
:type file: str
:param nvar: the number of frames of the movie file we want to plot
starting from the last line
:type nvar: int
:param lastvar: the index of the last timesteps to be read
:type lastvar: int
:param cm: the name of the color map
:type cm: str
:param levels: the number of contour levels
:type levels: int
:param cut: adjust the contour extrema to max(abs(data))*cut
:type cut: float
:param iplot: a boolean to specify if one wants to plot or not the
results
:type iplot: bool
:param avg: time average of the different forces
:type avg: bool
:param precision: precision of the input file, np.float32 for single
precision, np.float64 for double precision
:type precision: str
:param datadir: working directory
:type datadir: str
"""
if file is None:
dat = glob.glob('TO_mov.*')
str1 = 'Which TO movie do you want ?\n'
for k, movie in enumerate(dat):
str1 += ' {}) {}\n'.format(k+1, movie)
index = int(input(str1))
try:
filename = dat[index-1]
except IndexError:
print('Non valid index: {} has been chosen instead'.format(dat[0]))
filename = dat[0]
else:
filename = file
filename = os.path.join(datadir, filename)
# Get the version
self._get_version(filename)
# Read the movie file
infile = npfile(filename, endian='B')
# Header
self._read_header(infile, 0, precision)
if self.n_fields == 8:
self.l_phase = True
else:
self.l_phase = False
# Get the number of lines
self._get_nlines(datadir, filename, nvar, lastvar)
# Shape
self._get_data_shape(precision, allocate=False)
self.time = np.zeros(self.nvar, precision)
self.asVphi = np.zeros((self.nvar, self.n_theta_max, self.n_r_max),
precision)
self.rey = np.zeros_like(self.asVphi)
self.adv = np.zeros_like(self.asVphi)
self.visc = np.zeros_like(self.asVphi)
self.lorentz = np.zeros_like(self.asVphi)
self.coriolis = np.zeros_like(self.asVphi)
self.dtVp = np.zeros_like(self.asVphi)
if self.l_phase:
self.penalty = np.zeros_like(self.asVphi)
# Read the data
# If one skip the beginning, nevertheless read but do not store
for i in range(self.var2-self.nvar):
if self.version == 2:
n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
movieDipLon, movieDipStrength, movieDipStrengthGeo = \
infile.fort_read(precision)
else:
t_movieS = infile.fort_read(precision)
for ll in range(self.n_fields):
dat = infile.fort_read(precision, shape=self.shape, order='F')
for k in range(self.nvar):
if self.version == 2:
n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
movieDipLon, movieDipStrength, movieDipStrengthGeo \
= infile.fort_read(precision)
else:
t_movieS = infile.fort_read(precision)
self.time[k] = t_movieS
self.asVphi[k, ...] = infile.fort_read(precision, shape=self.shape,
order='F')
self.rey[k, ...] = infile.fort_read(precision, shape=self.shape,
order='F')
self.adv[k, ...] = infile.fort_read(precision, shape=self.shape,
order='F')
self.visc[k, ...] = infile.fort_read(precision, shape=self.shape,
order='F')
self.lorentz[k, ...] = infile.fort_read(precision, shape=self.shape,
order='F')
self.coriolis[k, ...] = infile.fort_read(precision, shape=self.shape,
order='F')
self.dtVp[k, ...] = infile.fort_read(precision, shape=self.shape,
order='F')
if self.l_phase:
self.penalty[k, ...] = infile.fort_read(precision, shape=self.shape,
order='F')
if iplot:
cmap = plt.get_cmap(cm)
self.plot(cut, levels, avg, cmap)
[docs] def __add__(self, new):
"""
Built-in function to sum two TO movies
.. note:: So far this function only works for two TO movies with the
same grid sizes. At some point, we might introduce grid
extrapolation to allow any summation.
"""
out = copy.deepcopy(new)
if new.time[0] == self.time[-1]:
out.time = np.concatenate((self.time, new.time[1:]), axis=0)
out.asVphi = np.concatenate((self.asVphi, new.asVphi[1:, ...]),
axis=0)
out.rey = np.concatenate((self.rey, new.rey[1:, ...]), axis=0)
out.adv = np.concatenate((self.adv, new.adv[1:, ...]), axis=0)
out.visc = np.concatenate((self.visc, new.visc[1:, ...]), axis=0)
out.lorentz = np.concatenate((self.lorentz, new.lorentz[1:, ...]),
axis=0)
out.coriolis = np.concatenate((self.coriolis, new.coriolis[1:, ...]),
axis=0)
out.dtVp = np.concatenate((self.dtVp, new.dtVp[1:, ...]), axis=0)
if self.l_phase:
out.penalty = np.concatenate((self.penalty, new.penalty[1:, ...]),
axis=0)
else:
out.time = np.concatenate((self.time, new.time), axis=0)
out.asVphi = np.concatenate((self.asVphi, new.asVphi), axis=0)
out.rey = np.concatenate((self.rey, new.rey), axis=0)
out.adv = np.concatenate((self.adv, new.adv), axis=0)
out.visc = np.concatenate((self.visc, new.visc), axis=0)
out.lorentz = np.concatenate((self.lorentz, new.lorentz), axis=0)
out.coriolis = np.concatenate((self.coriolis, new.coriolis),
axis=0)
out.dtVp = np.concatenate((self.dtVp, new.dtVp), axis=0)
if self.l_phase:
out.penalty = np.concatenate((self.penalty, new.penalty), axis=0)
return out
[docs] def plot(self, cut=0.8, levs=16, avg=True, cmap='RdYlBu_r'):
"""
Plotting function
:param cut: adjust the contour extrema to max(abs(data))*cut
:type cut: float
:param levs: number of contour levels
:type levs: int
:param avg: when set to True, quantities are time-averaged
:type avg: bool
:param cmap: name of the colormap
:type cmap: str
"""
th = np.linspace(np.pi/2., -np.pi/2., self.n_theta_max)
rr, tth = np.meshgrid(self.radius, th)
xx = rr * np.cos(tth)
yy = rr * np.sin(tth)
xxout = rr.max() * np.cos(th)
yyout = rr.max() * np.sin(th)
xxin = rr.min() * np.cos(th)
yyin = rr.min() * np.sin(th)
fig = plt.figure(figsize=(20, 5))
fig.subplots_adjust(top=0.99, right=0.99, bottom=0.01, left=0.01,
wspace=0.01)
if avg:
cor = self.coriolis.mean(axis=0)
asVp = self.asVphi.mean(axis=0)
rey = self.rey.mean(axis=0)
adv = self.adv.mean(axis=0)
lor = self.lorentz.mean(axis=0)
vis = self.visc.mean(axis=0)
if self.l_phase:
pen = self.penalty.mean(axis=0)
dt = self.dtVp[:-1].mean(axis=0)
vmin = - max(abs(cor.max()), abs(cor.min()))
vmin = cut * vmin
vmax = -vmin
cs = np.linspace(vmin, vmax, levs)
vmin = - max(abs(asVp.max()), abs(asVp.min()))
vmin = cut * vmin
vmax = -vmin
csVp = np.linspace(vmin, vmax, levs)
ax = fig.add_subplot(181)
ax.axis('off')
im = ax.contourf(xx, yy, asVp, csVp, cmap=cmap, extend='both')
ax.plot(xxout, yyout, 'k-', lw=1.5)
ax.plot(xxin, yyin, 'k-', lw=1.5)
ax.text(0.05, 0., 'uphi', fontsize=20, horizontalalignment='left',
verticalalignment='center')
ax = fig.add_subplot(182)
ax.axis('off')
im = ax.contourf(xx, yy, adv, cs, cmap=cmap, extend='both')
ax.plot(xxout, yyout, 'k-', lw=1.5)
ax.plot(xxin, yyin, 'k-', lw=1.5)
ax.text(0.05, 0., 'Adv', fontsize=20, horizontalalignment='left',
verticalalignment='center')
ax = fig.add_subplot(183)
ax.axis('off')
im = ax.contourf(xx, yy, rey, cs, cmap=cmap, extend='both')
ax.plot(xxout, yyout, 'k-', lw=1.5)
ax.plot(xxin, yyin, 'k-', lw=1.5)
ax.text(0.05, 0., 'Rey', fontsize=20, horizontalalignment='left',
verticalalignment='center')
ax = fig.add_subplot(184)
ax.axis('off')
im = ax.contourf(xx, yy, vis, cs, cmap=cmap, extend='both')
ax.plot(xxout, yyout, 'k-', lw=1.5)
ax.plot(xxin, yyin, 'k-', lw=1.5)
ax.text(0.05, 0., 'Visc', fontsize=20, horizontalalignment='left',
verticalalignment='center')
ax = fig.add_subplot(185)
ax.axis('off')
im = ax.contourf(xx, yy, lor, cs, cmap=cmap, extend='both')
ax.plot(xxout, yyout, 'k-', lw=1.5)
ax.plot(xxin, yyin, 'k-', lw=1.5)
ax.text(0.05, 0., 'Lo.', fontsize=20, horizontalalignment='left',
verticalalignment='center')
ax = fig.add_subplot(186)
ax.axis('off')
im = ax.contourf(xx, yy, cor, cs, cmap=cmap, extend='both')
ax.plot(xxout, yyout, 'k-', lw=1.5)
ax.plot(xxin, yyin, 'k-', lw=1.5)
ax.text(0.05, 0., 'Cor.', fontsize=20, horizontalalignment='left',
verticalalignment='center')
ax = fig.add_subplot(187)
ax.axis('off')
balance = adv+cor+vis+lor+rey
if self.l_phase:
balance += pen
im = ax.contourf(xx, yy, balance, cs, cmap=cmap, extend='both')
ax.plot(xxout, yyout, 'k-', lw=1.5)
ax.plot(xxin, yyin, 'k-', lw=1.5)
ax.text(0.05, 0., 'sum', fontsize=20, horizontalalignment='left',
verticalalignment='center')
ax = fig.add_subplot(188)
ax.axis('off')
im = ax.contourf(xx, yy, dt, cs, cmap=cmap, extend='both')
ax.plot(xxout, yyout, 'k-', lw=1.5)
ax.plot(xxin, yyin, 'k-', lw=1.5)
ax.text(0.01, 0., 'dvp/dt', fontsize=20,
horizontalalignment='left', verticalalignment='center')
else:
plt.ion()
vmin = - max(abs(self.coriolis.max()), abs(self.coriolis.min()))
vmin = cut * vmin
vmax = -vmin
cs = np.linspace(vmin, vmax, levs)
vmin = - max(abs(self.asVphi.max()), abs(self.asVphi.min()))
vmin = cut * vmin
vmax = -vmin
cs1 = np.linspace(vmin, vmax, levs)
for k in range(self.nvar-1): # avoid last dvp/dt which is wrong
bal = self.asVphi[k, ...]+self.adv[k, ...]+self.rey[k, ...] + \
self.visc[k, ...]+self.lorentz[k, ...] + \
self.coriolis[k, ...]
if k == 0:
ax1 = fig.add_subplot(181)
ax1.axis('off')
im = ax1.contourf(xx, yy, self.asVphi[k, ...], cs1,
cmap=cmap, extend='both')
ax1.plot(xxout, yyout, 'k-', lw=1.5)
ax1.plot(xxin, yyin, 'k-', lw=1.5)
ax1.text(0.05, 0., 'uphi', fontsize=20,
horizontalalignment='left',
verticalalignment='center')
ax2 = fig.add_subplot(182)
ax2.axis('off')
im = ax2.contourf(xx, yy, self.adv[k, ...], cs,
cmap=cmap, extend='both')
ax2.plot(xxout, yyout, 'k-', lw=1.5)
ax2.plot(xxin, yyin, 'k-', lw=1.5)
ax2.text(0.05, 0., 'Adv', fontsize=20,
horizontalalignment='left',
verticalalignment='center')
ax3 = fig.add_subplot(183)
ax3.axis('off')
im = ax3.contourf(xx, yy, self.rey[k, ...], cs,
cmap=cmap, extend='both')
ax3.plot(xxout, yyout, 'k-', lw=1.5)
ax3.plot(xxin, yyin, 'k-', lw=1.5)
ax3.text(0.05, 0., 'Rey', fontsize=20,
horizontalalignment='left',
verticalalignment='center')
ax4 = fig.add_subplot(184)
ax4.axis('off')
im = ax4.contourf(xx, yy, self.visc[k, ...], cs,
cmap=cmap, extend='both')
ax4.plot(xxout, yyout, 'k-', lw=1.5)
ax4.plot(xxin, yyin, 'k-', lw=1.5)
ax4.text(0.05, 0., 'Visc', fontsize=20,
horizontalalignment='left',
verticalalignment='center')
ax5 = fig.add_subplot(185)
ax5.axis('off')
im = ax5.contourf(xx, yy, self.lorentz[k, ...], cs,
cmap=cmap, extend='both')
ax5.plot(xxout, yyout, 'k-', lw=1.5)
ax5.plot(xxin, yyin, 'k-', lw=1.5)
ax5.text(0.05, 0., 'Lo.', fontsize=20,
horizontalalignment='left',
verticalalignment='center')
ax6 = fig.add_subplot(186)
ax6.axis('off')
im = ax6.contourf(xx, yy, self.coriolis[k, ...], cs,
cmap=cmap, extend='both')
ax6.plot(xxout, yyout, 'k-', lw=1.5)
ax6.plot(xxin, yyin, 'k-', lw=1.5)
ax6.text(0.05, 0., 'Cor.', fontsize=20,
horizontalalignment='left',
verticalalignment='center')
ax7 = fig.add_subplot(187)
ax7.axis('off')
im = ax7.contourf(xx, yy, bal, cs,
cmap=cmap, extend='both')
ax7.plot(xxout, yyout, 'k-', lw=1.5)
ax7.plot(xxin, yyin, 'k-', lw=1.5)
ax7.text(0.05, 0., 'sum', fontsize=20,
horizontalalignment='left',
verticalalignment='center')
ax8 = fig.add_subplot(188)
ax8.axis('off')
im = ax8.contourf(xx, yy, self.dtVp[k, ...], cs,
cmap=cmap, extend='both')
ax8.plot(xxout, yyout, 'k-', lw=1.5)
ax8.plot(xxin, yyin, 'k-', lw=1.5)
ax8.text(0.05, 0., 'dvp/dt', fontsize=20,
horizontalalignment='left',
verticalalignment='center')
man = plt.get_current_fig_manager()
man.canvas.draw()
else:
plt.cla()
ax1.contourf(xx, yy, self.asVphi[k, ...], cs1,
cmap=cmap, extend='both')
ax2.contourf(xx, yy, self.adv[k, ...], cs,
cmap=cmap, extend='both')
ax3.contourf(xx, yy, self.rey[k, ...], cs,
cmap=cmap, extend='both')
ax4.contourf(xx, yy, self.visc[k, ...], cs,
cmap=cmap, extend='both')
ax5.contourf(xx, yy, self.lorentz[k, ...], cs,
cmap=cmap, extend='both')
ax6.contourf(xx, yy, self.coriolis[k, ...], cs,
cmap=cmap, extend='both')
ax7.contourf(xx, yy, bal, cs,
cmap=cmap, extend='both')
ax8.contourf(xx, yy, self.dtVp[k, ...], cs,
cmap=cmap, extend='both')
plt.axis('off')
man.canvas.draw()
class MagicTOZ(MagicSetup):
"""
This class can be used to read the TOZ.TAG files produced by the TO outputs
>>> # read the content of TOZ_1.tag
>>> # where tag is the most recent file in the current directory
>>> toz = MagicTOZ(itoz=1)
>>> # read the content of TOZ_ave.test
>>> toz = MagicTOZ(tag='test', ave=True)
"""
def __init__(self, datadir='.', itoz=None, tag=None, precision=np.float32,
ave=False):
"""
:param datadir: current working directory
:type datadir: str
:param tag: file suffix (tag), if not specified the most recent one in
the current directory is chosen
:type tag: str
:param precision: single or double precision
:type precision: str
:param iplot: display the output plot when set to True (default is
True)
:type iplot: bool
:param ave: plot a time-averaged TOZ file when set to True
:type ave: bool
:param itoz: the number of the TOZ file you want to plot
:type itoz: int
"""
if ave:
self.name = 'TOZ_ave'
n_fields = 9
else:
self.name = 'TOZ_'
n_fields = 8
if tag is not None:
if itoz is not None:
file = '{}{}.{}'.format(self.name, itoz, tag)
filename = os.path.join(datadir, file)
else:
pattern = os.path.join(datadir, '{}*{}'.format(self.name, tag))
files = scanDir(pattern)
if len(files) != 0:
filename = files[-1]
else:
print('No such tag... try again')
return
if os.path.exists(os.path.join(datadir, 'log.{}'.format(tag))):
MagicSetup.__init__(self, datadir=datadir, quiet=True,
nml='log.{}'.format(tag))
else:
if itoz is not None:
pattern = os.path.join(datadir, '{}{}*'.format(self.name, itoz))
files = scanDir(pattern)
filename = files[-1]
else:
pattern = os.path.join(datadir, '{}*'.format(self.name))
files = scanDir(pattern)
filename = files[-1]
# Determine the setup
mask = re.compile(r'.*\.(.*)')
ending = mask.search(files[-1]).groups(0)[0]
if os.path.exists(os.path.join(datadir, 'log.{}'.format(ending))):
MagicSetup.__init__(self, datadir=datadir, quiet=True,
nml='log.{}'.format(ending))
infile = npfile(filename, endian='B')
if ave:
self.nSmax, self.omega_ic, self.omega_oc = infile.fort_read(precision)
else:
self.time, self.nSmax, self.omega_ic, \
self.omega_oc = infile.fort_read(precision)
self.nSmax = int(self.nSmax)
self.cylRad = infile.fort_read(precision)
self.zall = np.zeros((2*self.nSmax, self.nSmax), dtype=precision)
self.vp = np.zeros((2*self.nSmax, self.nSmax), dtype=precision)
self.dvp = np.zeros((2*self.nSmax, self.nSmax), dtype=precision)
self.Rstr = np.zeros((2*self.nSmax, self.nSmax), dtype=precision)
self.Astr = np.zeros((2*self.nSmax, self.nSmax), dtype=precision)
self.LF = np.zeros((2*self.nSmax, self.nSmax), dtype=precision)
self.Cor = np.zeros((2*self.nSmax, self.nSmax), dtype=precision)
self.visc = np.zeros((2*self.nSmax, self.nSmax), dtype=precision)
self.cylRadall = np.zeros_like(self.zall)
for i in range(2*self.nSmax):
self.cylRadall[i, :] = self.cylRad[:]
if n_fields == 9:
self.CL = np.zeros((2*self.nSmax, self.nSmax), dtype=precision)
for nS in range(self.nSmax):
nZmaxNS = int(infile.fort_read(precision))
print(nZmaxNS)
data = infile.fort_read(precision)
data = data.reshape((n_fields, nZmaxNS))
self.zall[:nZmaxNS, nS] = data[0, :]
self.vp[:nZmaxNS, nS] = data[1, :]
self.dvp[:nZmaxNS, nS] = data[2, :]
self.Rstr[:nZmaxNS, nS] = data[3, :]
self.Astr[:nZmaxNS, nS] = data[4, :]
self.LF[:nZmaxNS, nS] = data[5, :]
self.visc[:nZmaxNS, nS] = data[6, :]
self.Cor[:nZmaxNS, nS] = data[7, :]
if n_fields == 9:
self.CL[:nZmaxNS, nS] = data[8, :]
infile.close()
S = np.zeros_like(self.zall)
for i in range(2*self.nSmax):
S[i, :] = self.cylRad
[docs]class MagicTOHemi(MagicSetup):
"""
This class can be used to read and display z-integrated quantities
produced by the TO outputs. Those are basically the TO[s|n]hn.TAG files
>>> to = MagicTOHemi(hemi='n', iplot=True) # For the Northern hemisphere
"""
[docs] def __init__(self, datadir='.', hemi='n', tag=None, precision=np.float32,
iplot=False):
"""
:param datadir: current working directory
:type datadir: str
:param hemi: Northern or Southern hemisphere ('n' or 's')
:type hemi: str
:param tag: file suffix (tag), if not specified the most recent one in
the current directory is chosen
:type tag: str
:param precision: single or double precision
:type precision: str
:param iplot: display the output plot when set to True (default is
True)
:type iplot: bool
"""
if tag is not None:
pattern = os.path.join(datadir, 'TO{}hs.{}'.format(hemi, tag))
files = scanDir(pattern)
if len(files) != 0:
filename = files[-1]
else:
print('No such tag... try again')
return
if os.path.exists(os.path.join(datadir, 'log.{}'.format(tag))):
MagicSetup.__init__(self, datadir=datadir, quiet=True,
nml='log.{}'.format(tag))
else:
pattern = os.path.join(datadir, 'TO{}hs.*'.format(hemi))
files = scanDir(pattern)
filename = files[-1]
# Determine the setup
mask = re.compile(r'.*\.(.*)')
ending = mask.search(files[-1]).groups(0)[0]
if os.path.exists(os.path.join(datadir, 'log.{}'.format(ending))):
MagicSetup.__init__(self, datadir=datadir, quiet=True,
nml='log.{}'.format(ending))
if not os.path.exists(filename):
print('No such file')
return
infile = npfile(filename, endian='B')
self.nSmax = int(infile.fort_read(precision))
self.cylRad = infile.fort_read(precision)
self.height = np.zeros_like(self.cylRad)
ro = self.cylRad[0]
ri = ro - 1.
self.height[self.cylRad>=ri] = 2.*np.sqrt(ro**2-self.cylRad[self.cylRad>=ri]**2)
self.height[self.cylRad<ri] = np.sqrt(ro**2-self.cylRad[self.cylRad<ri]**2) - \
np.sqrt(ri**2-self.cylRad[self.cylRad<ri]**2)
ind = 0
while 1:
try:
data = infile.fort_read(precision)
if ind == 0:
self.time = np.r_[data[0]]
else:
self.time = np.append(self.time, data[0])
data = data[1:]
data = data.reshape((17, self.nSmax))
if ind == 0:
self.vp = data[0, :]
self.dvp = data[1, :]
self.ddvp = data[2, :]
self.vpr = data[3, :]
self.rstr = data[4, :]
self.astr = data[5, :]
self.LF = data[6, :]
self.viscstr = data[7, :]
self.tay = data[8, :]
self.Tau = data[9, :]
self.TauB = data[10, :]
self.BspdInt = data[11, :]
self.BpsdInt = data[12, :]
self.dTau = data[13, :]
self.dTTau = data[14, :]
self.dTTauB = data[15, :]
self.Bs2 = data[16, :]
else:
self.vp = np.vstack((self.vp, data[0, :]))
self.dvp = np.vstack((self.dvp, data[1, :]))
self.ddvp = np.vstack((self.ddvp, data[2, :]))
self.vpr = np.vstack((self.vpr, data[3, :]))
self.rstr = np.vstack((self.rstr, data[4, :]))
self.astr = np.vstack((self.astr, data[5, :]))
self.LF = np.vstack((self.LF, data[6, :]))
self.viscstr = np.vstack((self.viscstr, data[7, :]))
self.tay = np.vstack((self.tay, data[8, :]))
self.Tau = np.vstack((self.Tau, data[9, :]))
self.TauB = np.vstack((self.TauB, data[10, :]))
self.BspdInt = np.vstack((self.BspdInt, data[11, :]))
self.BpsdInt = np.vstack((self.BpsdInt, data[12, :]))
self.dTau = np.vstack((self.dTau, data[13, :]))
self.dTTau = np.vstack((self.dTTau, data[14, :]))
self.dTTauB = np.vstack((self.dTTauB, data[15, :]))
self.Bs2 = np.vstack((self.Bs2, data[16, :]))
ind += 1
except TypeError:
break
infile.close()
if iplot:
self.plot()
[docs] def __add__(self, new):
"""
This is an intrinsic '+' method to stack to TO[n|s]hs.TAG files
.. note:: So far, this only works if the grid resolution does not change
"""
out = copy.deepcopy(new)
if new.time[0] == self.time[-1]:
out.time = np.concatenate((self.time, new.time[1:]), axis=0)
out.vp = np.concatenate((self.vp, new.vp[1:, ...]), axis=0)
out.dvp = np.concatenate((self.dvp, new.dvp[1:, ...]), axis=0)
out.ddvp = np.concatenate((self.ddvp, new.ddvp[1:, ...]), axis=0)
out.tay = np.concatenate((self.tay, new.tay[1:, ...]), axis=0)
out.rstr = np.concatenate((self.rstr, new.rstr[1:, ...]), axis=0)
out.astr = np.concatenate((self.astr, new.astr[1:, ...]), axis=0)
out.viscstr = np.concatenate((self.viscstr, new.viscstr[1:, ...]), axis=0)
out.LF = np.concatenate((self.LF, new.LF[1:, ...]), axis=0)
out.Bs2 = np.concatenate((self.Bs2, new.Bs2[1:, ...]), axis=0)
else:
out.time = np.concatenate((self.time, new.time), axis=0)
out.vp = np.concatenate((self.vp, new.vp), axis=0)
out.dvp = np.concatenate((self.dvp, new.dvp), axis=0)
out.ddvp = np.concatenate((self.ddvp, new.ddvp), axis=0)
out.tay = np.concatenate((self.tay, new.tay), axis=0)
out.rstr = np.concatenate((self.rstr, new.rstr), axis=0)
out.astr = np.concatenate((self.astr, new.astr), axis=0)
out.viscstr = np.concatenate((self.viscstr, new.viscstr), axis=0)
out.LF = np.concatenate((self.LF, new.LF), axis=0)
out.Bs2 = np.concatenate((self.Bs2, new.Bs2), axis=0)
return out
[docs] def plot(self):
"""
Plotting function
"""
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(self.cylRad, self.vp.mean(axis=0))
ax.set_xlabel('Cylindrical radius')
ax.axhline(linestyle='--', linewidth=1.5, color='k')
ax.set_ylabel('Vphi')
#ax.set_xlim(self.cylRad[0], self.cylRad[-1])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(self.cylRad, self.rstr.mean(axis=0), label='Reynolds stress')
ax.plot(self.cylRad, self.astr.mean(axis=0),
label='Axi. Reynolds stress')
ax.plot(self.cylRad, self.LF.mean(axis=0), label='Lorentz force')
ax.plot(self.cylRad, self.viscstr.mean(axis=0), label='Viscous stress')
ax.set_xlabel('Cylindrical radius')
ax.set_ylabel('Integrated forces')
ax.legend(loc='upper left', frameon=False)
#ax.set_xlim(self.cylRad[0], self.cylRad[-1])
class MagicTaySphere(MagicSetup):
"""
This class can be used to read and display quantities
produced by the TO outputs. Those are basically the TaySphere.TAG files
>>> to = MagicTaySphere(iplot=True) # For the Northern hemisphere
>>> print(to.time, to.e_kin)
"""
def __init__(self, datadir='.', tag=None, precision=np.float32,
iplot=False):
"""
:param datadir: current working directory
:type datadir: str
:param tag: file suffix (tag), if not specified the most recent one in
the current directory is chosen
:type tag: str
:param precision: single or double precision
:type precision: str
:param iplot: display the output plot when set to True (default is
True)
:type iplot: bool
"""
if tag is not None:
pattern = os.path.join(datadir, 'TaySphere.{}'.format(tag))
files = scanDir(pattern)
if len(files) != 0:
filename = files[-1]
else:
print('No such tag... try again')
return
if os.path.exists(os.path.join(datadir, 'log.{}'.format(tag))):
MagicSetup.__init__(self, datadir=datadir, quiet=True,
nml='log.{}'.format(tag))
else:
pattern = os.path.join(datadir, 'TaySphere.*')
files = scanDir(pattern)
filename = files[-1]
# Determine the setup
mask = re.compile(r'.*\.(.*)')
ending = mask.search(files[-1]).groups(0)[0]
if os.path.exists(os.path.join(datadir, 'log.{}'.format(ending))):
MagicSetup.__init__(self, datadir=datadir, quiet=True,
nml='log.{}'.format(ending))
if not os.path.exists(filename):
print('No such file')
return
infile = npfile(filename, endian='B')
self.n_r_max = int(infile.fort_read(precision))
self.radius = infile.fort_read(precision)
ind = 0
while 1:
try:
data = infile.fort_read(precision)
if ind == 0:
self.time = np.r_[data[0]]
self.vpRMS = np.r_[data[1]]
self.vgRMS = np.r_[data[2]]
self.tayRMS = np.r_[data[3]]
self.taySRMS = np.r_[data[4]]
self.tayRRMS = np.r_[data[5]]
self.tayVRMS = np.r_[data[6]]
self.e_kin = np.r_[data[7]]
else:
self.time = np.append(self.time, data[0])
self.vpRMS = np.append(self.vgRMS, data[1])
self.vgRMS = np.append(self.vgRMS, data[2])
self.tayRMS = np.append(self.tayRMS, data[3])
self.taySRMS = np.append(self.taySRMS, data[4])
self.tayRRMS = np.append(self.tayRRMS, data[5])
self.tayVRMS = np.append(self.tayVRMS, data[6])
self.e_kin = np.append(self.e_kin, data[7])
data = data[8:]
data = data.reshape((8, self.n_r_max))
if ind == 0:
self.vp = data[0, :]
self.dvp = data[1, :]
self.rstr = data[2, :]
self.astr = data[3, :]
self.LF = data[4, :]
self.viscstr = data[5, :]
self.cor = data[6, :]
self.tay = data[7, :]
else:
self.vp = np.vstack((self.vp, data[0, :]))
self.dvp = np.vstack((self.dvp, data[1, :]))
self.rstr = np.vstack((self.rstr, data[2, :]))
self.astr = np.vstack((self.astr, data[3, :]))
self.LF = np.vstack((self.LF, data[4, :]))
self.viscstr = np.vstack((self.viscstr, data[5, :]))
self.cor = np.vstack((self.cor, data[6, :]))
self.tay = np.vstack((self.tay, data[7, :]))
ind += 1
except TypeError:
break
infile.close()
if iplot:
self.plot()
def plot(self):
"""
Plotting function
"""
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(self.radius, self.vp.mean(axis=0))
ax.set_xlabel('Radius')
ax.axhline(linestyle='--', linewidth=1.5, color='k')
ax.set_ylabel('Vphi')
ax.set_xlim(self.radius[0], self.radius[-1])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(self.radius, self.rstr.mean(axis=0), label='Reynolds stress')
ax.plot(self.radius, self.astr.mean(axis=0),
label='Axi. Reynolds stress')
ax.plot(self.radius, self.LF.mean(axis=0), label='Lorentz force')
ax.plot(self.radius, self.viscstr.mean(axis=0), label='Viscous stress')
ax.set_xlabel('Radius')
ax.set_ylabel('Forces')
ax.legend(loc='upper left', frameon=False)
ax.set_xlim(self.radius[0], self.radius[-1])
if __name__ == '__main__':
MagicTOZ(tag='start', itoz=1, ave=False, verbose=True)
MagicTOHemi(hemi='n', iplot=True)
file = 'TO_mov.test'
TOMovie(file=file)
plt.show()
