Skip to content
Snippets Groups Projects
Commit 7979b9a7 authored by Stefan Possanner's avatar Stefan Possanner
Browse files

Removed empty cells in notebook

parent 0069cafb
No related branches found
No related tags found
1 merge request!2Added License, delete old files, update README
Pipeline #154044 passed
Stage: test
Show whitespace changes
Inline Side-by-side
notebooks/plot_mhd_equil.ipynb
+ 4
− 39
View file @ 7979b9a7
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import os import os
import numpy as np import numpy as np
from gvec_to_python.reader.gvec_reader import create_GVEC_json from gvec_to_python.reader.gvec_reader import create_GVEC_json
from gvec_to_python import GVEC from gvec_to_python import GVEC
import gvec_to_python as _ import gvec_to_python as _
from matplotlib import pyplot as plt from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d from mpl_toolkits.mplot3d import axes3d
# give absolute paths to the files # give absolute paths to the files
pkg_path = _.__path__[0] pkg_path = _.__path__[0]
print('package path:', pkg_path) print('package path:', pkg_path)
# dat_file_in = os.path.join( # dat_file_in = os.path.join(
# pkg_path, 'testcases/ellipstell/GVEC_ellipStell_profile_update_State_0000_00010000.dat') # pkg_path, 'testcases/ellipstell/GVEC_ellipStell_profile_update_State_0000_00010000.dat')
# json_file_out = os.path.join( # json_file_out = os.path.join(
# pkg_path, 'testcases/ellipstell/GVEC_ellipStell_profile_update_State_0000_00010000.json') # pkg_path, 'testcases/ellipstell/GVEC_ellipStell_profile_update_State_0000_00010000.json')
dat_file_in = os.path.join( dat_file_in = os.path.join(
pkg_path, 'testcases/ellipstell_v2/newBC_E1D6_M6N6/GVEC_ELLIPSTELL_V2_State_0000_00200000.dat') pkg_path, 'testcases/ellipstell_v2/newBC_E1D6_M6N6/GVEC_ELLIPSTELL_V2_State_0000_00200000.dat')
json_file_out = os.path.join( json_file_out = os.path.join(
pkg_path, 'testcases/ellipstell_v2/newBC_E1D6_M6N6/GVEC_ELLIPSTELL_V2_State_0000_00200000.json') pkg_path, 'testcases/ellipstell_v2/newBC_E1D6_M6N6/GVEC_ELLIPSTELL_V2_State_0000_00200000.json')
create_GVEC_json(dat_file_in, json_file_out, with_spl_coef=False) create_GVEC_json(dat_file_in, json_file_out, with_spl_coef=False)
# main object (one without, the other one with pyccel kernels) # main object (one without, the other one with pyccel kernels)
gvec = GVEC(json_file_out, use_nfp=True) gvec = GVEC(json_file_out, use_nfp=False)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Plot profiles # Plot profiles
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
s = np.linspace(0, 1, 100) s = np.linspace(0, 1, 100)
fig = plt.figure(figsize=(13, 10)) fig = plt.figure(figsize=(13, 10))
# toroidal flux # toroidal flux
ax = fig.add_subplot(2, 2, 1) ax = fig.add_subplot(2, 2, 1)
ax.plot(s, gvec.profiles.profile(s, name='phi')) ax.plot(s, gvec.profiles.profile(s, name='phi'))
ax.set_xlabel('s') ax.set_xlabel('s')
ax.set_ylabel('$\Phi(s)$') ax.set_ylabel('$\Phi(s)$')
ax.set_title('Toroidal flux') ax.set_title('Toroidal flux')
# poloidal flux # poloidal flux
ax = fig.add_subplot(2, 2, 2) ax = fig.add_subplot(2, 2, 2)
ax.plot(s, gvec.profiles.profile(s, name='chi')) ax.plot(s, gvec.profiles.profile(s, name='chi'))
ax.set_xlabel('s') ax.set_xlabel('s')
ax.set_ylabel('$\chi(s)$') ax.set_ylabel('$\chi(s)$')
ax.set_title('Poloidal flux') ax.set_title('Poloidal flux')
# iota # iota
ax = fig.add_subplot(2, 2, 3) ax = fig.add_subplot(2, 2, 3)
ax.plot(s, gvec.profiles.profile(s, name='iota')) ax.plot(s, gvec.profiles.profile(s, name='iota'))
ax.set_xlabel('s') ax.set_xlabel('s')
ax.set_ylabel('$\iota(s)$') ax.set_ylabel('$\iota(s)$')
ax.set_title('Iota') ax.set_title('Iota')
# pressure # pressure
ax = fig.add_subplot(2, 2, 4) ax = fig.add_subplot(2, 2, 4)
ax.plot(s, gvec.profiles.profile(s, name='pressure')) ax.plot(s, gvec.profiles.profile(s, name='pressure'))
ax.set_xlabel('s') ax.set_xlabel('s')
ax.set_ylabel('$p(s)$') ax.set_ylabel('$p(s)$')
ax.set_title('Pressure') ax.set_title('Pressure')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Poloidal geometry # Poloidal geometry
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
gvec.mapping = 'unit' # use the mapping setter gvec.mapping = 'unit' # use the mapping setter
n_planes = 7 n_planes = 4
s = np.linspace(0, 1, 20) # radial coordinate in [0, 1] s = np.linspace(0, 1, 20) # radial coordinate in [0, 1]
u = np.linspace(0, 1, 49) # poloidal angle in [0, 1] u = np.linspace(0, 1, 49) # poloidal angle in [0, 1]
v = np.linspace(0, 1, n_planes) # toroidal angle in [0, 1] v = np.linspace(0, 1, n_planes) # toroidal angle in [0, 1]
x, y, z = gvec.f(s, u, v) x, y, z = gvec.f(s, u, v)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(13, np.ceil(n_planes/2) * 6.5)) fig = plt.figure(figsize=(13, np.ceil(n_planes/2) * 6.5))
# poloidal plane grid # poloidal plane grid
for n in range(v.size): for n in range(v.size):
xp = x[:, :, n].squeeze() xp = x[:, :, n].squeeze()
yp = y[:, :, n].squeeze() yp = y[:, :, n].squeeze()
zp = z[:, :, n].squeeze() zp = z[:, :, n].squeeze()
rp = np.sqrt(xp**2 + yp**2) rp = np.sqrt(xp**2 + yp**2)
ax = fig.add_subplot(int(np.ceil(n_planes/2)), 2, n + 1) ax = fig.add_subplot(int(np.ceil(n_planes/2)), 2, n + 1)
for i in range(rp.shape[0]): for i in range(rp.shape[0]):
for j in range(rp.shape[1] - 1): for j in range(rp.shape[1] - 1):
if i < rp.shape[0] - 1: if i < rp.shape[0] - 1:
ax.plot([rp[i, j], rp[i + 1, j]], [zp[i, j], zp[i + 1, j]], 'b', linewidth=.6) ax.plot([rp[i, j], rp[i + 1, j]], [zp[i, j], zp[i + 1, j]], 'b', linewidth=.6)
if j < rp.shape[1] - 1: if j < rp.shape[1] - 1:
ax.plot([rp[i, j], rp[i, j + 1]], [zp[i, j], zp[i, j + 1]], 'b', linewidth=.6) ax.plot([rp[i, j], rp[i, j + 1]], [zp[i, j], zp[i, j + 1]], 'b', linewidth=.6)
ax.set_xlabel('r') ax.set_xlabel('r')
ax.set_ylabel('z') ax.set_ylabel('z')
ax.set_title('Poloidal plane at $\\theta$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp)) ax.set_title('Poloidal plane at $\\theta$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp))
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Jacobian determinant # Jacobian determinant
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
det_df = gvec.det_df(s, u, v) det_df = gvec.det_df(s, u, v)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(13, np.ceil(n_planes/2) * 6.5)) fig = plt.figure(figsize=(13, np.ceil(n_planes/2) * 6.5))
# poloidal plane grid # poloidal plane grid
for n in range(v.size): for n in range(v.size):
xp = x[:, :, n].squeeze() xp = x[:, :, n].squeeze()
yp = y[:, :, n].squeeze() yp = y[:, :, n].squeeze()
zp = z[:, :, n].squeeze() zp = z[:, :, n].squeeze()
rp = np.sqrt(xp**2 + yp**2) rp = np.sqrt(xp**2 + yp**2)
detp = det_df[:, :, n].squeeze() detp = det_df[:, :, n].squeeze()
ax = fig.add_subplot(int(np.ceil(n_planes/2)), 2, n + 1) ax = fig.add_subplot(int(np.ceil(n_planes/2)), 2, n + 1)
map = ax.contourf(rp, zp, detp, 30) map = ax.contourf(rp, zp, detp, 30)
ax.set_xlabel('r') ax.set_xlabel('r')
ax.set_ylabel('z') ax.set_ylabel('z')
ax.set_title('Jacobian determinant at $\\theta$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp)) ax.set_title('Jacobian determinant at $\\theta$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp))
fig.colorbar(map, ax=ax, location='right') fig.colorbar(map, ax=ax, location='right')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Top view # Top view
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
s = np.linspace(0, 1, 20) # radial coordinate in [0, 1] s = np.linspace(0, 1, 20) # radial coordinate in [0, 1]
u = np.linspace(0, 1, 3) # poloidal angle in [0, 1] u = np.linspace(0, 1, 3) # poloidal angle in [0, 1]
v = np.linspace(0, 1, 69) # toroidal angle in [0, 1] v = np.linspace(0, 1, 69) # toroidal angle in [0, 1]
x, y, z = gvec.f(s, u, v) x, y, z = gvec.f(s, u, v)
fig = plt.figure(figsize=(13, 2 * 6.5)) fig = plt.figure(figsize=(13, 2 * 6.5))
ax = fig.add_subplot() ax = fig.add_subplot()
# poloidal plane grid # poloidal plane grid
for m in range(2): for m in range(2):
xp = x[:, m, :].squeeze() xp = x[:, m, :].squeeze()
yp = y[:, m, :].squeeze() yp = y[:, m, :].squeeze()
zp = z[:, m, :].squeeze() zp = z[:, m, :].squeeze()
for i in range(xp.shape[0]): for i in range(xp.shape[0]):
for j in range(xp.shape[1] - 1): for j in range(xp.shape[1] - 1):
if i < xp.shape[0] - 1: if i < xp.shape[0] - 1:
ax.plot([xp[i, j], xp[i + 1, j]], [yp[i, j], yp[i + 1, j]], 'b', linewidth=.6) ax.plot([xp[i, j], xp[i + 1, j]], [yp[i, j], yp[i + 1, j]], 'b', linewidth=.6)
if j < xp.shape[1] - 1: if j < xp.shape[1] - 1:
if i == 0: if i == 0:
ax.plot([xp[i, j], xp[i, j + 1]], [yp[i, j], yp[i, j + 1]], 'r', linewidth=1) ax.plot([xp[i, j], xp[i, j + 1]], [yp[i, j], yp[i, j + 1]], 'r', linewidth=1)
else: else:
ax.plot([xp[i, j], xp[i, j + 1]], [yp[i, j], yp[i, j + 1]], 'b', linewidth=.6) ax.plot([xp[i, j], xp[i, j + 1]], [yp[i, j], yp[i, j + 1]], 'b', linewidth=.6)
ax.set_xlabel('x') ax.set_xlabel('x')
ax.set_ylabel('y') ax.set_ylabel('y')
ax.axis('equal') ax.axis('equal')
ax.set_title('Stellarator top view') ax.set_title('Stellarator top view')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Pressure # Pressure
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
n_planes = 4 n_planes = 4
s = np.linspace(0, 1, 20) # radial coordinate in [0, 1] s = np.linspace(0, 1, 20) # radial coordinate in [0, 1]
u = np.linspace(0, 1, 49) # poloidal angle in [0, 1] u = np.linspace(0, 1, 49) # poloidal angle in [0, 1]
v = np.linspace(0, 1, n_planes) # toroidal angle in [0, 1] v = np.linspace(0, 1, n_planes) # toroidal angle in [0, 1]
p, xyz = gvec.p_cart(s, u, v) p, xyz = gvec.p_cart(s, u, v)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(15, np.ceil(n_planes/2) * 6.5)) fig = plt.figure(figsize=(15, np.ceil(n_planes/2) * 6.5))
# poloidal plane grid # poloidal plane grid
for n in range(v.size): for n in range(v.size):
xp = xyz[0][:, :, n].squeeze() xp = xyz[0][:, :, n].squeeze()
yp = xyz[1][:, :, n].squeeze() yp = xyz[1][:, :, n].squeeze()
zp = xyz[2][:, :, n].squeeze() zp = xyz[2][:, :, n].squeeze()
rp = np.sqrt(xp**2 + yp**2) rp = np.sqrt(xp**2 + yp**2)
pp = p[:, :, n].squeeze() pp = p[:, :, n].squeeze()
ax = fig.add_subplot(int(np.ceil(n_planes/2)), 2, n + 1) ax = fig.add_subplot(int(np.ceil(n_planes/2)), 2, n + 1)
map = ax.contourf(rp, zp, pp, 30) map = ax.contourf(rp, zp, pp, 30)
ax.set_xlabel('r') ax.set_xlabel('r')
ax.set_ylabel('z') ax.set_ylabel('z')
ax.set_title('Pressure at $\\theta$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp)) ax.set_title('Pressure at $\\theta$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp))
fig.colorbar(map, ax=ax, location='right') fig.colorbar(map, ax=ax, location='right')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Magnetic field # Magnetic field
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
n_planes = 4 n_planes = 4
s = np.linspace(0.0001, 1, 20) # radial coordinate in [0, 1] s = np.linspace(0.0001, 1, 20) # radial coordinate in [0, 1]
u = np.linspace(0, 1, 39) # poloidal angle in [0, 1] u = np.linspace(0, 1, 39) # poloidal angle in [0, 1]
v = np.linspace(0, 1, n_planes) # toroidal angle in [0, 1] v = np.linspace(0, 1, n_planes) # toroidal angle in [0, 1]
b, xyz = gvec.b_cart(s, u, v) b, xyz = gvec.b_cart(s, u, v)
df = gvec.df(s,u,v) df = gvec.df(s,u,v)
abs_b = np.sqrt(b[0]**2 + b[1]**2 + b[2]**2) abs_b = np.sqrt(b[0]**2 + b[1]**2 + b[2]**2)
print(abs_b.shape) print(abs_b.shape)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(15, np.ceil(n_planes/2) * 6.5)) fig = plt.figure(figsize=(15, np.ceil(n_planes/2) * 6.5))
for n in range(v.size): for n in range(v.size):
xp = xyz[0][:, :, n].squeeze() xp = xyz[0][:, :, n].squeeze()
yp = xyz[1][:, :, n].squeeze() yp = xyz[1][:, :, n].squeeze()
zp = xyz[2][:, :, n].squeeze() zp = xyz[2][:, :, n].squeeze()
rp = np.sqrt(xp**2 + yp**2) rp = np.sqrt(xp**2 + yp**2)
ab = abs_b[:, :, n].squeeze() ab = abs_b[:, :, n].squeeze()
ax = fig.add_subplot(int(np.ceil(n_planes/2)), 2, n + 1) ax = fig.add_subplot(int(np.ceil(n_planes/2)), 2, n + 1)
map = ax.contourf(rp, zp, ab, 30) map = ax.contourf(rp, zp, ab, 30)
ax.set_xlabel('r') ax.set_xlabel('r')
ax.set_ylabel('z') ax.set_ylabel('z')
ax.set_title('Magnetic field strength at $\\theta$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp)) ax.set_title('Magnetic field strength at $\\theta$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp))
fig.colorbar(map, ax=ax, location='right') fig.colorbar(map, ax=ax, location='right')
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(13, np.ceil(n_planes/2) * 6.5)) fig = plt.figure(figsize=(13, np.ceil(n_planes/2) * 6.5))
for n in range(v.size): for n in range(v.size):
xp = xyz[0][:, :, n].squeeze() xp = xyz[0][:, :, n].squeeze()
yp = xyz[1][:, :, n].squeeze() yp = xyz[1][:, :, n].squeeze()
zp = xyz[2][:, :, n].squeeze() zp = xyz[2][:, :, n].squeeze()
rp = np.sqrt(xp**2 + yp**2) rp = np.sqrt(xp**2 + yp**2)
# plot theta derivative of mapping, # plot theta derivative of mapping,
dxdu = df[:, :, n,0,1].squeeze() dxdu = df[:, :, n,0,1].squeeze()
dydu = df[:, :, n,1,1].squeeze() dydu = df[:, :, n,1,1].squeeze()
dzdu = df[:, :, n,2,1].squeeze() dzdu = df[:, :, n,2,1].squeeze()
phi = 2*np.pi*v[n] / gvec.domain.nfp phi = 2*np.pi*v[n] / gvec.domain.nfp
drdu = dxdu*np.cos(phi) - dydu*np.sin(phi) drdu = dxdu*np.cos(phi) - dydu*np.sin(phi)
ax = fig.add_subplot(int(np.ceil(n_planes/2)), 2, n + 1) ax = fig.add_subplot(int(np.ceil(n_planes/2)), 2, n + 1)
ax.quiver(rp, zp, drdu, dzdu, scale=40) ax.quiver(rp, zp, drdu, dzdu, scale=40)
ax.set_xlabel('r') ax.set_xlabel('r')
ax.set_ylabel('z') ax.set_ylabel('z')
ax.set_title('Theta derivative of mapping at $\\phi$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp)) ax.set_title('Theta derivative of mapping at $\\phi$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp))
#fig.colorbar(map, ax=ax, location='right') #fig.colorbar(map, ax=ax, location='right')
``` ```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment