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

Update notebook

parent c0911911
No related branches found
No related tags found
1 merge request!5Added current density j
%% 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)
#gvec_file='testcases/ellipstell/newBC_E1D6_M6N6/GVEC_ELLIPSTELL_E1D6_M6N6_State_0000_00200000' #gvec_file='testcases/ellipstell/newBC_E1D6_M6N6/GVEC_ELLIPSTELL_E1D6_M6N6_State_0000_00200000'
gvec_file='testcases/ellipstell/newBC_E4D6_M6N6/GVEC_ELLIPSTELL_E4D6_M6N6_State_0001_00200000' gvec_file='testcases/ellipstell/newBC_E4D6_M6N6/GVEC_ELLIPSTELL_E4D6_M6N6_State_0001_00200000'
dat_file_in = os.path.join(pkg_path, gvec_file+'.dat') dat_file_in = os.path.join(pkg_path, gvec_file+'.dat')
json_file_out = os.path.join(pkg_path, gvec_file+'.json') json_file_out = os.path.join(pkg_path, gvec_file+'.json')
create_GVEC_json(dat_file_in, json_file_out) create_GVEC_json(dat_file_in, json_file_out)
# 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=False) 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 = 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]
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
gvec.mapping = 'unit' # use the mapping setter
n_planes = 4
s = np.linspace(0, 1, 20) # radial coordinate 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]
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
gvec.mapping = 'unit' # use the mapping setter
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)
```
%% Cell type:code id: tags:
``` python
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
gvec.mapping = 'unit' # use the mapping setter
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
gvec.mapping = 'unit' # use the mapping setter
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 $\\phi$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp)) ax.set_title('Magnetic field strength 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: %% 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:markdown id: tags: %% Cell type:markdown id: tags:
# Current density # Current density
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
gvec.mapping = 'gvec' gvec.mapping = 'gvec'
n_planes = 4 n_planes = 4
s = np.linspace(0.2, 1, 10) # radial coordinate in [0, 1] s = np.linspace(0.2, 1, 10) # radial coordinate in [0, 1]
u = np.linspace(0, 1, 19) # poloidal angle in [0, 1] u = np.linspace(0, 2*np.pi, 19) # poloidal angle in [0, 1]
v = np.linspace(0, 1, n_planes) # toroidal angle in [0, 1] v = np.linspace(0, 2*np.pi, n_planes) # toroidal angle in [0, 1]
j, xyz = gvec.j_cart(s, u, v) j, xyz = gvec.j_cart(s, u, v)
abs_j = np.sqrt(j[0]**2 + j[1]**2 + j[2]**2) abs_j = np.sqrt(j[0]**2 + j[1]**2 + j[2]**2)
print(abs_j.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_j[:, :, n].squeeze() ab = abs_j[:, :, 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('Current (absolute value) at $\\phi$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp)) ax.set_title('Current (absolute value) 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: %% Cell type:code id: tags:
``` python ``` python
gvec.assert_equil(s, u, v) gvec.assert_equil(s, u, v)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
p = gvec.profiles.profile(s, u, v, name='pressure', der=None)
p_s = gvec.profiles.profile(s, u, v, name='pressure', der='s')
j = gvec.j2(s, u, v)
b = gvec.bv(s, u, v)
abs_j = np.sqrt(j[0]**2 + j[1]**2 + j[2]**2)
abs_b = np.sqrt(b[0]**2 + b[1]**2 + b[2]**2)
```
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(15, np.ceil(n_planes/2) * 6.5))
ss, uu, vv = gvec.domain.prepare_args(s, u, v, sparse=False)
print(ss.shape, uu.shape, vv.shape)
comp = 1
for n in range(v.size):
sp = ss[:, :, n].squeeze()
up = uu[:, :, n].squeeze()
vp = vv[:, :, n].squeeze()
ji = j[comp][:, :, n].squeeze()
ax = fig.add_subplot(int(np.ceil(n_planes/2)), 2, n + 1)
map = ax.contourf(sp, up, ji, 30)
ax.set_xlabel('s')
ax.set_ylabel('u')
ax.set_title('Current component {1} at $\\phi$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp, comp))
fig.colorbar(map, ax=ax, location='right')
```
%% Cell type:code id: tags:
``` python
``` ```
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment