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

Notebook update

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