Skip to content
Snippets Groups Projects
Commit 4eed92f7 authored by Florian Hindenlang's avatar Florian Hindenlang Committed by Stefan Possanner
Browse files

fix in grad btheta, fixes current computation

parent a4046513
Branches
Tags
1 merge request!10fix in grad btheta, fixes current computation
Show whitespace changes
Inline Side-by-side
notebooks/plot_mhd_equil.ipynb
+ 88
− 16
View file @ 4eed92f7
%% 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,unit_tor_domain="one-fp",use_pyccel=True) gvec = GVEC(json_file_out,unit_tor_domain="one-fp",use_pyccel=True)
``` ```
%% 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
def get_uvfac(mapping): def get_uvfac(mapping):
'''Depending on the mapping, gives back the factor for u in [0,1] (either 1 or 2pi) '''Depending on the mapping, gives back the factor for u in [0,1] (either 1 or 2pi)
and v in[0,1] scaled to always represent one field period''' and v in[0,1] scaled to always represent one field period'''
if(mapping=='unit' or mapping=='unit_pest'): if(mapping=='unit' or mapping=='unit_pest'):
ufac=1. ufac=1.
# 'unit' mapping can include a tor_fraction. We want to always visualize one field period. # 'unit' mapping can include a tor_fraction. We want to always visualize one field period.
vfac = 1./(gvec.domain.tor_fraction*gvec.nfp) vfac = 1./(gvec.domain.tor_fraction*gvec.nfp)
else: #'gvec',"wo_hmap",'pest' else: #'gvec',"wo_hmap",'pest'
ufac=2*np.pi ufac=2*np.pi
vfac=2*np.pi/gvec.nfp vfac=2*np.pi/gvec.nfp
return ufac,vfac return ufac,vfac
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
gvec.mapping = 'pest' # use the mapping setter gvec.mapping = 'gvec' # 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]
ufac,vfac=get_uvfac(gvec.mapping) # compute factors to always visualize one field period ufac,vfac=get_uvfac(gvec.mapping) # compute factors to always visualize one field period
a1 = np.linspace(0, 1*ufac, 39) # poloidal angle in [0, 1] a1 = np.linspace(0, 1*ufac, 39) # poloidal angle in [0, 1]
a2 = np.linspace(0, 1*vfac, n_planes,endpoint=False) a2 = np.linspace(0, 1*vfac, n_planes,endpoint=False)
x, y, z = gvec.f(s, a1, a2) x, y, z = gvec.f(s, a1, a2)
``` ```
%% 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(a2.size):
if(gvec.mapping=="wo_hmap"): if(gvec.mapping=="wo_hmap"):
rp = x[:, :, n].squeeze() rp = x[:, :, n].squeeze()
zp = y[:, :, n].squeeze() zp = y[:, :, n].squeeze()
else: else:
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.axis('equal') ax.axis('equal')
ax.set_title('Poloidal plane at $\\zeta$={0:4.3f} $\pi$'.format(a2[n]/vfac*2/gvec.nfp)) ax.set_title('Poloidal plane at $\\zeta$={0:4.3f} $\pi$'.format(a2[n]/vfac*2/gvec.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 = 'gvec' # use the mapping setter gvec.mapping = 'gvec' # 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]
ufac,vfac=get_uvfac(gvec.mapping) # compute factors to always visualize one field period ufac,vfac=get_uvfac(gvec.mapping) # compute factors to always visualize one field period
u = np.linspace(0, 1*ufac, 49) # poloidal angle u = np.linspace(0, 1*ufac, 49) # poloidal angle
v = np.linspace(0, 1*vfac, n_planes,endpoint=False) #toroidal angle v = np.linspace(0, 1*vfac, n_planes,endpoint=False) #toroidal angle
x, y, z = gvec.f(s, u, v) x, y, z = gvec.f(s, u, v)
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.axis('equal') ax.axis('equal')
ax.set_title('Jacobian determinant at $\\zeta$={0:4.3f} $\pi$'.format(v[n]/vfac*2/gvec.nfp)) ax.set_title('Jacobian determinant at $\\zeta$={0:4.3f} $\pi$'.format(v[n]/vfac*2/gvec.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 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: %% Cell type:code id: tags:
``` python ``` 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 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]
ufac,vfac=get_uvfac(gvec.mapping) # compute factors to always visualize one field period ufac,vfac=get_uvfac(gvec.mapping) # compute factors to always visualize one field period
u = np.linspace(0, 1*ufac, 49) # poloidal angle u = np.linspace(0, 1*ufac, 49) # poloidal angle
v = np.linspace(0, 1*vfac, n_planes,endpoint=False) #toroidal angle v = np.linspace(0, 1*vfac, n_planes,endpoint=False) #toroidal angle
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.axis('equal') ax.axis('equal')
ax.set_title('Pressure at $\\zeta$={0:4.3f} $\pi$'.format(v[n]/vfac * 2 / gvec.nfp)) ax.set_title('Pressure at $\\zeta$={0:4.3f} $\pi$'.format(v[n]/vfac * 2 / gvec.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 = 'gvec' # use the mapping setter gvec.mapping = 'gvec' # use the mapping setter
n_planes = 4 n_planes = 4
s = np.linspace(0.0001, 1, 13) # radial coordinate in [0, 1] s = np.linspace(0.0001, 1, 13) # radial coordinate in [0, 1]
ufac,vfac=get_uvfac(gvec.mapping) # compute factors to always visualize one field period ufac,vfac=get_uvfac(gvec.mapping) # compute factors to always visualize one field period
u = np.linspace(0, 1*ufac, 49) # poloidal angle u = np.linspace(0, 1*ufac, 49) # poloidal angle
v = np.linspace(0, 1*vfac, n_planes,endpoint=False) #toroidal angle v = np.linspace(0, 1*vfac, n_planes,endpoint=False) #toroidal angle
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)
if(gvec.mapping=="wo_hmap"): if(gvec.mapping=="wo_hmap"):
abs_b = np.sqrt(b[0]**2 + b[1]**2 + xyz[0]**2*b[2]**2) abs_b = np.sqrt(b[0]**2 + b[1]**2 + xyz[0]**2*b[2]**2)
else: else:
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):
if(gvec.mapping=="wo_hmap"): if(gvec.mapping=="wo_hmap"):
rp = xyz[0][:, :, n].squeeze() rp = xyz[0][:, :, n].squeeze()
zp = xyz[1][:, :, n].squeeze() zp = xyz[1][:, :, n].squeeze()
else: else:
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.axis('equal') ax.axis('equal')
ax.set_title('Magnetic field strength at $\\zeta$={0:4.3f} $\pi$'.format(v[n]/vfac * 2 / gvec.nfp)) ax.set_title('Magnetic field strength at $\\zeta$={0:4.3f} $\pi$'.format(v[n]/vfac * 2 / gvec.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()
# 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()
zeta = 2*np.pi/ gvec.nfp *v[n]/vfac zeta = 2*np.pi/ gvec.nfp *v[n]/vfac
drdu = dxdu*np.cos(zeta) - dydu*np.sin(zeta) drdu = dxdu*np.cos(zeta) - dydu*np.sin(zeta)
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=20) ax.quiver(rp, zp, drdu, dzdu, scale=20)
ax.set_xlabel('r') ax.set_xlabel('R')
ax.set_ylabel('z') ax.set_ylabel('Z')
ax.axis('equal') ax.axis('equal')
ax.set_title('Theta derivative of mapping at $\\zeta$={0:4.3f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp)) ax.set_title('Theta derivative of mapping at $\\zeta$={0:4.3f} $\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
def bcart(s,th,ze):
b,x=gvec.b_cart(s,th,ze)
return b
def jcart_FD(s,th,ze):
eps=1.0e-8
b_cart = bcart(s,th,ze)
b_cart_s_eps = bcart(s+eps,th,ze)
b_cart_th_eps = bcart(s,th+eps,ze)
b_cart_ze_eps = bcart(s,th,ze+eps)
gradb=np.zeros((3,3))
for d in [0,1,2]:
b_ds =(b_cart_s_eps[d] -b_cart[d])/eps
b_dth=(b_cart_th_eps[d]-b_cart[d])/eps
b_dze=(b_cart_ze_eps[d]-b_cart[d])/eps
gradb[d] = gvec.df_inv(s,th,ze).T@np.array([b_ds,b_dth,b_dze])
j_cart_FD=np.zeros(3)
for d in [0,1,2]:
j_cart_FD[d] = gradb[(d+2)%3][(d+1)%3]-gradb[(d+1)%3][(d+2)%3]
return j_cart_FD
gvec.mapping = 'gvec'
s,th,ze = [0.49,0.5,0.3]
mu0=4.0e-7*np.pi
j_ex ,xyz = gvec.j_cart(s, th, ze)
j_FD = jcart_FD(s,th,ze)/mu0
print('j_FD ',j_FD)
print('j_ex',np.array(j_ex))
print('j_ex-j_FD',np.array(j_ex)-j_FD)
assert(np.allclose(j_ex,j_FD,rtol=1e-6,atol=1e-6/mu0))
```
%% Cell type:code id: tags:
``` python
gvec.mapping = 'gvec' gvec.mapping = 'gvec'
n_planes = 4 n_planes = 4
s = np.linspace(0.001, 1, 10) # radial coordinate in [0, 1] s = np.linspace(0.001, 1, 10) # radial coordinate in [0, 1]
ufac,vfac=get_uvfac(gvec.mapping) # compute factors to always visualize one field period ufac,vfac=get_uvfac(gvec.mapping) # compute factors to always visualize one field period
u = np.linspace(0, 1*ufac, 49) # poloidal angle u = np.linspace(0, 1*ufac, 49) # poloidal angle
v = np.linspace(0, 1*vfac, n_planes,endpoint=False) #toroidal angle v = np.linspace(0, 1*vfac, n_planes,endpoint=False) #toroidal angle
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)
``` ```
%% 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.axis('equal') ax.axis('equal')
ax.set_title('Current (absolute value) at $\\zeta$={0:4.3f} $\pi$'.format(v[n]/vfac * 2 / gvec.nfp)) ax.set_title('Current (absolute value) at $\\zeta$={0:4.3f} $\pi$'.format(v[n]/vfac * 2 / gvec.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 = gvec.profiles.profile(s, u, v, name='pressure', der=None)
p_s = gvec.profiles.profile(s, u, v, name='pressure', der='s') p_s = gvec.profiles.profile(s, u, v, name='pressure', der='s')
j = gvec.j2(s, u, v) j = gvec.j2(s, u, v)
b = gvec.bv(s, u, v) b = gvec.bv(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)
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)
``` ```
%% 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))
ss, uu, vv = gvec.domain.prepare_args(s, u, v, sparse=False) ss, uu, vv = gvec.domain.prepare_args(s, u, v, sparse=False)
print(ss.shape, uu.shape, vv.shape) print(ss.shape, uu.shape, vv.shape)
comp = 1 comp = 0
for n in range(v.size): for n in range(v.size):
sp = ss[:, :, n].squeeze() sp = ss[:, :, n].squeeze()
up = uu[:, :, n].squeeze() up = uu[:, :, n].squeeze()
vp = vv[:, :, n].squeeze() vp = vv[:, :, n].squeeze()
ji = j[comp][:, :, n].squeeze() ji = j[comp][:, :, 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(sp, up, ji, 30) map = ax.contourf(sp, up, ji, 30)
ax.set_xlabel('s') ax.set_xlabel('s')
ax.set_ylabel('u') ax.set_ylabel('u')
ax.axis('equal')
ax.set_title('Current component {1} at $\\phi$={0:4.1f} $\pi$'.format(v[n] * 2 / gvec.domain.nfp, comp)) 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') fig.colorbar(map, ax=ax, location='right')
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` 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 = 0
for n in range(v.size):
sp = ss[:, :, n].squeeze()
up = uu[:, :, n].squeeze()
vp = vv[:, :, n].squeeze()
if(comp==0):
fbe = (p_s - (j[1] * b[2] - j[2] * b[1]))[:, :, n].squeeze()
elif(comp==1):
fbe = (j[2] * b[0] - j[0] * b[2])[:, :, n].squeeze()
elif(comp==2):
fbe = (j[0] * b[1] - j[1] * b[0])[:, :, n].squeeze()
ax = fig.add_subplot(int(np.ceil(n_planes/2)), 2, n + 1)
map = ax.contourf(sp, up, fbe, 30)
ax.set_xlabel('s')
ax.set_ylabel('u')
ax.set_title('force balance error in comp {1} at $\\zeta$={0:4.1f} $\pi$'.format(v[n]/vfac * 2 / gvec.nfp, comp))
fig.colorbar(map, ax=ax, location='right')
```
%% Cell type:code id: tags:
``` python
``` ```
......
src/gvec_to_python/GVEC_functions.py
+ 2
− 2
View file @ 4eed92f7
...@@ -581,11 +581,11 @@ class GVEC: ...@@ -581,11 +581,11 @@ class GVEC:
term2 = - self._b_small_th(s, th, ze) * self.domain.dJ_gvec(s, th, ze, der='s') / det_df term2 = - self._b_small_th(s, th, ze) * self.domain.dJ_gvec(s, th, ze, der='s') / det_df
db_ds = (term1 + term2) / det_df db_ds = (term1 + term2) / det_df
term1 = dchi - dphi * self._lambda_thze(s, th, ze) term1 = - dphi * self._lambda_thze(s, th, ze)
term2 = - self._b_small_th(s, th, ze) * self.domain.dJ_gvec(s, th, ze, der='th') / det_df term2 = - self._b_small_th(s, th, ze) * self.domain.dJ_gvec(s, th, ze, der='th') / det_df
db_dt = (term1 + term2) / det_df db_dt = (term1 + term2) / det_df
term1 = dchi - dphi * self._lambda_zeze(s, th, ze) term1 = - dphi * self._lambda_zeze(s, th, ze)
term2 = - self._b_small_th(s, th, ze) * self.domain.dJ_gvec(s, th, ze, der='ze') / det_df term2 = - self._b_small_th(s, th, ze) * self.domain.dJ_gvec(s, th, ze, der='ze') / det_df
db_dz = (term1 + term2) / det_df db_dz = (term1 + term2) / det_df
......
src/gvec_to_python/geometry/domain.py
+ 15
− 17
View file @ 4eed92f7
...@@ -880,12 +880,14 @@ class GVEC_domain: ...@@ -880,12 +880,14 @@ class GVEC_domain:
for n in range(3): for n in range(3):
# meshgrid evaluation # meshgrid evaluation
if dX.ndim == 5: if dX.ndim == 5:
tmp += [ddX[0][:, :, :, 0, n] * dX2_dt + ddX[1][:, :, :, 1, n] * dX1_ds tmp += [ddX[n][:, :, :, 0, 0] * dX2_dt + ddX[n][:, :, :, 1, 1] * dX1_ds
- ddX[0][:, :, :, 1, n] * dX1_dt - ddX[1][:, :, :, 0, n] * dX2_ds] - ddX[n][:, :, :, 1, 0] * dX1_dt - ddX[n][:, :, :, 0, 1] * dX2_ds]
# single point evaluation # single point evaluation
# X1ds/di * dX2_dt + dX1dt/di * dX1_ds
# -X2ds/di * dX1_dt - dX2dt/di * dX2_ds
else: else:
tmp += [ddX[0][0, n] * dX2_dt + ddX[1][1, n] * dX1_ds tmp += [ddX[n][0, 0] * dX2_dt + ddX[n][1, 1] * dX1_ds
- ddX[0][1, n] * dX1_dt - ddX[1][0, n] * dX2_ds] - ddX[n][1, 0] * dX1_dt - ddX[n][0, 1] * dX2_ds]
return tmp[0], tmp[1], tmp[2] return tmp[0], tmp[1], tmp[2]
...@@ -1119,10 +1121,6 @@ class GVEC_domain: ...@@ -1119,10 +1121,6 @@ class GVEC_domain:
s.ndim, v.ndim) s.ndim, v.ndim)
if s.ndim == 1: if s.ndim == 1:
s, u, v = np.meshgrid(s, u, v, indexing='ij', sparse=sparse) s, u, v = np.meshgrid(s, u, v, indexing='ij', sparse=sparse)
elif s.ndim == 3:
is_meshgrid=(np.all(np.minimum.reduce(s,axis=(1,2))==np.maximum.reduce(s,axis=(1,2))) and
np.all(np.minimum.reduce(u,axis=(0,2))==np.maximum.reduce(u,axis=(0,2))) and
np.all(np.minimum.reduce(v,axis=(0,1))==np.maximum.reduce(v,axis=(0,1))))
return s, u, v return s, u, v
......
tests/test_gvec_equil.py
+ 79
− 0
View file @ 4eed92f7
...@@ -212,6 +212,85 @@ def test_values(use_pyccel,gvec_file,unit_tor_domain): ...@@ -212,6 +212,85 @@ def test_values(use_pyccel,gvec_file,unit_tor_domain):
assert np.allclose(a_unit_pest[1], a_gvec[1]) assert np.allclose(a_unit_pest[1], a_gvec[1])
assert np.allclose(a_unit_pest[2], a_gvec[2]) assert np.allclose(a_unit_pest[2], a_gvec[2])
@pytest.mark.parametrize('use_pyccel,gvec_file', [(True ,"testcases/ellipstell/newBC_E4D6_M6N6/GVEC_ELLIPSTELL_E4D6_M6N6_State_0001_00200000")])
@pytest.mark.parametrize('mapping', ['gvec', 'pest', 'unit', 'unit_pest'])
def test_derivs_with_FD(mapping,use_pyccel,gvec_file):
'''Compare exact derivatives of functions to their finite difference counterpart.'''
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 _
def do_FD(func,s,th,ze,der=None):
''' compute s,th,ze derivative of input function func'''
eps=1.0e-8
f0 = func(s,th,ze)
df_ds = (func(s+eps,th,ze) - f0 )/ eps
df_dth = (func(s,th+eps,ze) - f0 )/ eps
df_dze = (func(s,th,ze+eps) - f0 )/ eps
if(der=="s"):
return df_ds
elif(der=="th"):
return df_dth
elif(der=="ze"):
return df_dze
else:
return df_ds,df_dth,df_dze
def bth(s,th,ze):
''' computes B^theta=b^theta/det_df'''
return gvec._b_small_th(s, th, ze) / gvec.domain.det_df_gvec(s, th, ze)
def bze(s,th,ze):
''' computes B^zeta=b^zeta/det_df'''
return gvec._b_small_ze(s, th, ze) / gvec.domain.det_df_gvec(s, th, ze)
# give absolute paths to the files
pkg_path = _.__path__[0]
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
gvec = GVEC(json_file_out, mapping=mapping, use_pyccel=use_pyccel)
print('gvec.mapping: ', gvec.mapping)
# test single-point float evaluation
s = .45
if 'unit' not in mapping:
th = 0.77*2*np.pi
ze = 0.33*2*np.pi
else:
th = .77
ze = .33
#check derivative of Jacobian of gvec mapping
dJf_FD = do_FD(gvec.domain.det_df_gvec,s,th,ze)
dJf_ex = [gvec.domain.dJ_gvec(s,th,ze,der=deriv) for deriv in ["s","th","ze"]]
assert(np.allclose(dJf_FD,dJf_ex,rtol=1e-6,atol=1e-4))
# check derivative of metric tensor g of the gvec mapping
for deriv in ["s","th","ze"]:
dg_FD = do_FD(gvec.domain.g_gvec,s,th,ze,der=deriv)
dg_ex = gvec.domain.dg_gvec(s,th,ze,der=deriv)
assert(np.allclose(dg_FD,dg_ex,rtol=1e-6,atol=1e-4))
# check derivative of B^theta
dbth_FD = do_FD(bth,s,th,ze)
dbth_ex = gvec._grad_b_theta(s, th, ze)
assert(np.allclose(dbth_FD,dbth_ex,rtol=1e-6,atol=1e-4))
# check derivative of B^zeta
dbze_FD = do_FD(bze,s,th,ze)
dbze_ex = gvec._grad_b_zeta(s, th, ze)
#print("diff dbze_FD-dbze_ex",dbze_FD-np.array(dbze_ex))
assert(np.allclose(dbze_FD,dbze_ex,rtol=1e-6,atol=1e-4))
if __name__ == '__main__': if __name__ == '__main__':
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment