scripts/gray/visual.py: degrade instead of failing on missing files
This commit is contained in:
Michele Guerini Rocco
parent
9442356797
commit
0ee7b683fd
GPG Key ID:
BFBAF4C975F76450
@ -39,6 +39,20 @@ def cli_args() -> argparse.Namespace:
|
||||
return parser.parse_args()
|
||||
|
||||
|
||||
def placeholder(ax: plt.Axes, text: str):
|
||||
'''
|
||||
Puts a placeholder text at the center of an axis
|
||||
'''
|
||||
left, width = .25, .5
|
||||
bottom, height = .25, .5
|
||||
right = left + width
|
||||
top = bottom + height
|
||||
|
||||
ax.text(0.5*(left+right), 0.5*(top+bottom), text,
|
||||
horizontalalignment='center',
|
||||
verticalalignment='center')
|
||||
|
||||
|
||||
def align_yaxis(*axes: [plt.Axes]):
|
||||
'''
|
||||
Aligns the origins of two axes
|
||||
@ -116,30 +130,33 @@ def plot_poloidal(inputs: Path, outputs: Path, ax: plt.Axes):
|
||||
ax.set_ylabel('$z$ / m')
|
||||
|
||||
# load flux surfaces
|
||||
surfaces = gray.read_table(outputs / 'flux-surfaces.71.txt')
|
||||
surfaces = surfaces.reshape(-1, int(surfaces['i'].max()))
|
||||
try:
|
||||
surfaces = gray.read_table(outputs / 'flux-surfaces.71.txt')
|
||||
surfaces = surfaces.reshape(-1, int(surfaces['i'].max()))
|
||||
|
||||
# plot plasma boundary
|
||||
bound = np.argmin(1 - surfaces['ψ_n'][:, 0])
|
||||
ax.plot(surfaces[bound]['R'], surfaces[bound]['z'],
|
||||
c='xkcd:slate grey', label='plasma boundary')
|
||||
# plot plasma boundary
|
||||
bound = np.argmin(1 - surfaces['ψ_n'][:, 0])
|
||||
ax.plot(surfaces[bound]['R'], surfaces[bound]['z'],
|
||||
c='xkcd:slate grey', label='plasma boundary')
|
||||
|
||||
# manual contourplot of ψ(R,z)
|
||||
for i, surf in enumerate(surfaces[:bound]):
|
||||
ax.plot(surf['R'], surf['z'], c='xkcd:grey',
|
||||
label='flux surfaces' if i == 0 else '')
|
||||
# manual contourplot of ψ(R,z)
|
||||
for i, surf in enumerate(surfaces[:bound]):
|
||||
ax.plot(surf['R'], surf['z'], c='xkcd:grey',
|
||||
label='flux surfaces' if i == 0 else '')
|
||||
|
||||
# label the rational surfaces
|
||||
labels = ['q=2', 'q=3/2', 'q=1']
|
||||
for label, surf in zip(labels, surfaces[1+bound:][::-1]):
|
||||
line = np.vstack([surf['R'], surf['z']]).T
|
||||
cont = ContourSet(ax, [1], [[line]], filled=False,
|
||||
colors='xkcd:ocean blue',
|
||||
linestyles='--')
|
||||
ax.clabel(cont, [1], inline=True,
|
||||
fontsize=10, fmt={1: label})
|
||||
ax.plot(np.nan, np.nan, '--', color='xkcd:ocean blue',
|
||||
label='rational surfaces')
|
||||
# label the rational surfaces
|
||||
labels = ['q=2', 'q=3/2', 'q=1']
|
||||
for label, surf in zip(labels, surfaces[1+bound:][::-1]):
|
||||
line = np.vstack([surf['R'], surf['z']]).T
|
||||
cont = ContourSet(ax, [1], [[line]], filled=False,
|
||||
colors='xkcd:ocean blue',
|
||||
linestyles='--')
|
||||
ax.clabel(cont, [1], inline=True,
|
||||
fontsize=10, fmt={1: label})
|
||||
ax.plot(np.nan, np.nan, '--', color='xkcd:ocean blue',
|
||||
label='rational surfaces')
|
||||
except FileNotFoundError:
|
||||
pass
|
||||
|
||||
# load limiter
|
||||
limiter = gray.get_limiter(conf, inputs)
|
||||
@ -186,12 +203,15 @@ def plot_toroidal(inputs: Path, outputs: Path, ax: plt.Axes):
|
||||
limiter = gray.get_limiter(conf, inputs)
|
||||
|
||||
# plot plasma boundary
|
||||
surfaces = gray.read_table(outputs / 'flux-surfaces.71.txt')
|
||||
boundary = surfaces[np.isclose(surfaces['ψ_n'], 1, 1e-3)]
|
||||
try:
|
||||
surfaces = gray.read_table(outputs / 'flux-surfaces.71.txt')
|
||||
boundary = surfaces[np.isclose(surfaces['ψ_n'], 1, 1e-3)]
|
||||
|
||||
# plot plasma boundary
|
||||
plot_circle(ax, radius=boundary['R'].min(), color='xkcd:slate gray')
|
||||
plot_circle(ax, radius=boundary['R'].max(), color='xkcd:slate gray')
|
||||
# plot plasma boundary
|
||||
plot_circle(ax, radius=boundary['R'].min(), color='xkcd:slate gray')
|
||||
plot_circle(ax, radius=boundary['R'].max(), color='xkcd:slate gray')
|
||||
except FileNotFoundError:
|
||||
pass
|
||||
|
||||
# plot limiter
|
||||
if limiter[0].size > 0:
|
||||
@ -374,36 +394,42 @@ def plot_inputs(inputs: Path, outputs: Path, axes: [plt.Axes]):
|
||||
'''
|
||||
from matplotlib.colors import TwoSlopeNorm
|
||||
|
||||
maps = gray.read_table(outputs / 'inputs-maps.72.txt')
|
||||
conf = gray.read_conf(inputs / 'gray.ini')
|
||||
|
||||
# filter valid points
|
||||
maps = maps[maps['ψ_n'] > 0]
|
||||
flux = maps['R'], maps['z'], maps['ψ_n']
|
||||
try:
|
||||
maps = gray.read_table(outputs / 'inputs-maps.72.txt')
|
||||
# filter valid points
|
||||
maps = maps[maps['ψ_n'] > 0]
|
||||
flux = maps['R'], maps['z'], maps['ψ_n']
|
||||
except FileNotFoundError:
|
||||
maps = None
|
||||
placeholder(axes['B'], 'inputs-maps.72.txt not found')
|
||||
|
||||
# contour levels
|
||||
ψ_max = maps['ψ_n'].max() if maps is not None else 1.2
|
||||
inner_levels = np.linspace(0, 0.99, 8)
|
||||
outer_levels = np.linspace(0.9991, maps['ψ_n'].max()*0.99, 6)
|
||||
outer_levels = np.linspace(0.9991, ψ_max*0.99, 6)
|
||||
all_levels = np.concatenate([inner_levels, outer_levels])
|
||||
|
||||
# contour style
|
||||
norm = TwoSlopeNorm(vmin=0, vcenter=1, vmax=maps['ψ_n'].max())
|
||||
norm = TwoSlopeNorm(vmin=0, vcenter=1, vmax=ψ_max)
|
||||
cmap = plt.cm.RdYlGn_r
|
||||
borders = dict(colors='xkcd:grey', linestyles='-', linewidths=0.8)
|
||||
|
||||
# interpolated equilibrium
|
||||
interp = axes['B']
|
||||
interp.set_title('poloidal flux', loc='right')
|
||||
interp.set_xlabel('$R$ / m')
|
||||
interp.set_ylabel('$z$ / m')
|
||||
interp.tricontourf(*flux, levels=all_levels, norm=norm, cmap=cmap)
|
||||
interp.tricontour(*flux, levels=all_levels, **borders)
|
||||
interp.tricontour(*flux, levels=[1], colors='xkcd:slate gray')
|
||||
interp.plot(np.nan, np.nan, c='xkcd:slate gray', label='plasma boundary')
|
||||
if maps is not None:
|
||||
interp = axes['B']
|
||||
interp.set_title('poloidal flux', loc='right')
|
||||
interp.set_xlabel('$R$ / m')
|
||||
interp.set_ylabel('$z$ / m')
|
||||
interp.tricontourf(*flux, levels=all_levels, norm=norm, cmap=cmap)
|
||||
interp.tricontour(*flux, levels=all_levels, **borders)
|
||||
interp.tricontour(*flux, levels=[1], colors='xkcd:slate gray')
|
||||
interp.plot(np.nan, np.nan, c='xkcd:slate gray', label='plasma boundary')
|
||||
|
||||
# add limiter
|
||||
conf = gray.read_conf(inputs / 'gray.ini')
|
||||
limiter = gray.get_limiter(conf, inputs)
|
||||
interp.plot(*limiter, c='xkcd:black', label='limiter')
|
||||
# add limiter
|
||||
limiter = gray.get_limiter(conf, inputs)
|
||||
interp.plot(*limiter, c='xkcd:black', label='limiter')
|
||||
|
||||
# original MHD equilibrium
|
||||
orig = axes['A']
|
||||
@ -418,14 +444,18 @@ def plot_inputs(inputs: Path, outputs: Path, axes: [plt.Axes]):
|
||||
orig.plot(*eqdsk.limiter, c='xkcd:black')
|
||||
orig.plot(*eqdsk.boundary, c='xkcd:slate gray')
|
||||
else:
|
||||
orig.axis('off')
|
||||
placeholder(orig, 'not an EQDSK')
|
||||
|
||||
# colorbar
|
||||
bar = plt.colorbar(plt.cm.ScalarMappable(norm, cmap), cax=axes['C'])
|
||||
bar.set_label(label='normalised poloidal flux', labelpad=10)
|
||||
|
||||
# Plasma radial profiles
|
||||
profiles = gray.read_table(outputs / 'kinetic-profiles.55.txt')
|
||||
try:
|
||||
profiles = gray.read_table(outputs / 'kinetic-profiles.55.txt')
|
||||
except FileNotFoundError:
|
||||
placeholder(axes['D'], 'kinetic-profiles.55.txt not found')
|
||||
return
|
||||
|
||||
axes['D'].set_title('plasma profiles', loc='right')
|
||||
add_alt_axis(axes['D'], profiles['ρ_t'], profiles['ψ_n']**0.5,
|
||||
|
||||
Loading…
Reference in New Issue
Block a user