Particle tracing#

In this tutorial you will build and compare particle-tracing simulations with the toy models Vlasov and GuidingCenter.

By the end, you will know how to:

  1. change the geometry,

  2. choose different particle loading strategies,

  3. add a static background magnetic field.

Part 1: Particles in a cylinder#

As in Tutorial 2, we recreate the setup directly in the notebook so each option is visible and easy to modify. If you prefer starting from default launch files in a terminal, run:

struphy params Vlasov
struphy params GuidingCenter

Here we set the simulation domain \(\Omega\) to a cylinder and compare two sampling strategies for particle positions:

  • uniform in logical space \([0,1]^3 = F^{-1}(\Omega)\),

  • uniform directly in the physical cylinder \(\Omega\).

Start by importing the API and the model:

[1]:

from struphy import BaseUnits, EnvironmentOptions, Time
from struphy import equils
from struphy import domains
from struphy import grids
from struphy import DerhamOptions
from struphy import maxwellians
from struphy import (BoundaryParameters,
                     LoadingParameters,
                     WeightsParameters,
                     SortingParameters,
                     SavingParameters,
                     )
from struphy import Simulation

# import model
from struphy.models import Vlasov

Step 1: define separate output folders#

We will run two comparable simulations and save them in different output folders. Set the folder names through environment variables:

[2]:

# light-weight model instances
model = Vlasov()
model_2 = Vlasov()

# environment options
env = EnvironmentOptions()
env_2 = EnvironmentOptions(sim_folder="sim_2")

Step 2: set shared simulation options#

Keep the generic options identical in both runs so the loading strategy is the main difference. We perform only one time step, use a cylindrical geometry, and leave the equilibrium empty:

[3]:

# time stepping
time_opts = Time(dt=0.2, Tend=0.2)

# geometry
a1 = 1e-6 # to avoid polar singularity, we cut a small hole around the axis.
a2 = 5.0
Lz = 20.0
domain = domains.HollowCylinder(a1=a1, a2=a2, Lz=Lz)

# fluid equilibrium (can be used as part of initial conditions)
equil = None

# grid
grid = grids.TensorProductGrid()

# derham options
derham_opts = DerhamOptions()

Step 3: instantiate the first simulation#

Create the first simulation object with the baseline particle loading:

[4]:

sim = Simulation(model,
                 env=env,
                 time_opts=time_opts,
                 domain=domain,
                 equil=equil,
                 grid=grid,
                 derham_opts=derham_opts,)

Struphy provides a convenient way to clone a simulation and change only selected parameters.

Use that pattern here to create a second simulation with modified model/environment options:

[5]:

sim_2 = sim.spawn_sister(model= model_2, env=env_2)

Step 4: inspect the domain#

Before changing the loading options, verify the simulation domain visually:

[6]:

domain.show()
../../_images/_collections_tutorials_tutorial_particle_tracing_11_0.png

For the second simulation, set spatial="disc" in the loading parameters.

This samples particles uniformly on the cylinder cross section, which makes it easier to compare against logical-space loading:

[7]:

# species parameters
loading_params = LoadingParameters(Np=1000)
loading_params_2 = LoadingParameters(Np=1000, spatial="disc")

weights_params = WeightsParameters()
boundary_params = BoundaryParameters()
saving_params = SavingParameters(n_markers=1.0)

model.kinetic_ions.set_markers(
    loading_params=loading_params,
    weights_params=weights_params,
    boundary_params=boundary_params,
    saving_params=saving_params
)
model_2.kinetic_ions.set_markers(
    loading_params=loading_params_2,
    weights_params=weights_params,
    boundary_params=boundary_params,
    saving_params=saving_params
)

Use identical propagator options and initial conditions in both runs so that only the loading choice changes the result:

[8]:

# propagator options
model.propagators.push_vxb.options = model.propagators.push_vxb.Options()
model.propagators.push_eta.options = model.propagators.push_eta.Options()

model_2.propagators.push_vxb.options = model_2.propagators.push_vxb.Options()
model_2.propagators.push_eta.options = model_2.propagators.push_eta.Options()

[9]:

# initial conditions (background + perturbation)
perturbation = None
background = maxwellians.Maxwellian3D(n=(1.0, perturbation))

model.kinetic_ions.var.add_background(background)
model_2.kinetic_ions.var.add_background(background)

Step 5: run simulation A#

Execute the first simulation:

[10]:

out = sim.run()

WARNING: Class "BasisProjectionOperators" called with degree=(1, 1, 1) (interpolation of piece-wise constants should be avoided).

Time stepping: 100%|██████████| 1/1 [00:00<00:00, 151.55step/s]

Step 6: run simulation B#

Now execute the second simulation with disc-based spatial loading:

[11]:

out_2 = sim_2.run()

WARNING: Class "BasisProjectionOperators" called with degree=(1, 1, 1) (interpolation of piece-wise constants should be avoided).

Time stepping: 100%|██████████| 1/1 [00:00<00:00, 154.83step/s]

Step 7: compare initial particle distributions#

Load the generated data and plot initial particle positions on a cylinder cross section.

This plot is the key check for understanding how loading in logical space differs from loading in physical space:

[12]:

from matplotlib import pyplot as plt

fig = plt.figure(figsize=(10, 6))

# out.evaluate("<species>/orbits") returns an xarray.Dataset with one (t, marker) variable per
# saved quantity, e.g. orbits.x, orbits.v1, orbits.weight; print(orbits) lists them all
orbits = out.evaluate("kinetic_ions/orbits")
orbits_uni = out_2.evaluate("kinetic_ions/orbits")

# initial marker positions (first saved time step)
pos = orbits.isel(t=0)
pos_uni = orbits_uni.isel(t=0)

plt.subplot(1, 2, 1)
plt.scatter(pos.x, pos.y, s=2.0)
circle1 = plt.Circle((0, 0), a2, color="k", fill=False)
ax = plt.gca()
ax.add_patch(circle1)
ax.set_aspect("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.title("sim_1: draw uniform in logical space")

plt.subplot(1, 2, 2)
plt.scatter(pos_uni.x, pos_uni.y, s=2.0)
circle2 = plt.Circle((0, 0), a2, color="k", fill=False)
ax = plt.gca()
ax.add_patch(circle2)
ax.set_aspect("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.title("sim_2: draw uniform on disc")


No post-processed data in /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_1, processing with default options (call out.pproc(...) to choose them)

Post-processing path /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_1

No feec fields found in hdf5 file, skipping post-processing of fields.
Evaluation of 1000 marker orbits for kinetic_ions

100%|██████████| 2/2 [00:00<00:00, 227.69it/s]


No post-processed data in /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_2, processing with default options (call out.pproc(...) to choose them)

Post-processing path /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_2

No feec fields found in hdf5 file, skipping post-processing of fields.
Evaluation of 1000 marker orbits for kinetic_ions

100%|██████████| 2/2 [00:00<00:00, 224.79it/s]

[12]:

Text(0.5, 1.0, 'sim_2: draw uniform on disc')
../../_images/_collections_tutorials_tutorial_particle_tracing_22_5.png

Optional: quasi-uniform marker loading (Sobol)#

For a more homogeneous particle representation with limited particle count, use Sobol-based loading.

Set loading="sobol_standard" or loading="sobol_antithetic" in LoadingParameters:

[13]:

env_3 = EnvironmentOptions(sim_folder="sim_3")
model_3 = Vlasov()
sim_3 = Simulation.spawn_sister(sim, env=env_3, model=model_3)

# species parameters
loading_params_3 = LoadingParameters(Np=1000, spatial="disc", loading="sobol_standard")

model_3.kinetic_ions.set_markers(
    loading_params=loading_params_3,
    weights_params=weights_params,
    boundary_params=boundary_params,
    saving_params=saving_params,
)

# propagator options
model_3.propagators.push_vxb.options = model_3.propagators.push_vxb.Options()
model_3.propagators.push_eta.options = model_3.propagators.push_eta.Options()

# initial conditions (background + perturbation)
model_3.kinetic_ions.var.add_background(background)

out_3 = sim_3.run()

fig = plt.figure(figsize=(15, 6))

orbits_standard = out_3.evaluate("kinetic_ions/orbits")

plt.subplot(1, 3, 1)
plt.scatter(orbits.x[0], orbits.y[0], s=2.0)
circle1 = plt.Circle((0, 0), a2, color="k", fill=False)
ax = plt.gca()
ax.add_patch(circle1)
ax.set_aspect("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.title("sim_1: draw uniform in logical space")

plt.subplot(1, 3, 2)
plt.scatter(orbits_uni.x[0], orbits_uni.y[0], s=2.0)
circle2 = plt.Circle((0, 0), a2, color="k", fill=False)
ax = plt.gca()
ax.add_patch(circle2)
ax.set_aspect("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.title("sim_2: draw uniform on disc")

plt.subplot(1, 3, 3)
plt.scatter(orbits_standard.x[0], orbits_standard.y[0], s=2.0)
circle3 = plt.Circle((0, 0), a2, color="k", fill=False)
ax = plt.gca()
ax.add_patch(circle3)
ax.set_aspect("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.title("sim_3: draw using sobol_standard")

WARNING: Class "BasisProjectionOperators" called with degree=(1, 1, 1) (interpolation of piece-wise constants should be avoided).

Time stepping: 100%|██████████| 1/1 [00:00<00:00, 127.93step/s]


No post-processed data in /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_3, processing with default options (call out.pproc(...) to choose them)

Post-processing path /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_3

No feec fields found in hdf5 file, skipping post-processing of fields.
Evaluation of 1000 marker orbits for kinetic_ions


100%|██████████| 2/2 [00:00<00:00, 212.43it/s]

[13]:

Text(0.5, 1.0, 'sim_3: draw using sobol_standard')
../../_images/_collections_tutorials_tutorial_particle_tracing_24_5.png

Optional: antithetic Sobol loading#

Repeat the same setup with loading="sobol_antithetic" if you want a variance-reduced counterpart to standard Sobol loading:

[14]:

env_3 = EnvironmentOptions(sim_folder="sim_3")
model_3 = Vlasov()
sim_3 = Simulation.spawn_sister(sim, env=env_3, model=model_3)

# species parameters
loading_params_3 = LoadingParameters(Np=1000, spatial="disc", loading="sobol_standard")

model_3.kinetic_ions.set_markers(
    loading_params=loading_params_3,
    weights_params=weights_params,
    boundary_params=boundary_params,
    saving_params=saving_params,
)

# propagator options
model_3.propagators.push_vxb.options = model_3.propagators.push_vxb.Options()
model_3.propagators.push_eta.options = model_3.propagators.push_eta.Options()

# initial conditions (background + perturbation)
model_3.kinetic_ions.var.add_background(background)

out_3 = sim_3.run()

fig = plt.figure(figsize=(15, 6))

orbits_standard = out_3.evaluate("kinetic_ions/orbits")

plt.subplot(1, 3, 1)
plt.scatter(orbits.x[0], orbits.y[0], s=2.0)
circle1 = plt.Circle((0, 0), a2, color="k", fill=False)
ax = plt.gca()
ax.add_patch(circle1)
ax.set_aspect("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.title("sim_1: draw uniform in logical space")

plt.subplot(1, 3, 2)
plt.scatter(orbits_uni.x[0], orbits_uni.y[0], s=2.0)
circle2 = plt.Circle((0, 0), a2, color="k", fill=False)
ax = plt.gca()
ax.add_patch(circle2)
ax.set_aspect("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.title("sim_2: draw uniform on disc")

plt.subplot(1, 3, 3)
plt.scatter(orbits_standard.x[0], orbits_standard.y[0], s=2.0)
circle3 = plt.Circle((0, 0), a2, color="k", fill=False)
ax = plt.gca()
ax.add_patch(circle3)
ax.set_aspect("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.title("sim_3: draw using sobol_standard")

WARNING: Class "BasisProjectionOperators" called with degree=(1, 1, 1) (interpolation of piece-wise constants should be avoided).

Time stepping: 100%|██████████| 1/1 [00:00<00:00, 160.73step/s]


No post-processed data in /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_3, processing with default options (call out.pproc(...) to choose them)

Post-processing path /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_3

No feec fields found in hdf5 file, skipping post-processing of fields.
Evaluation of 1000 marker orbits for kinetic_ions


100%|██████████| 2/2 [00:00<00:00, 222.78it/s]

[14]:

Text(0.5, 1.0, 'sim_3: draw using sobol_standard')
../../_images/_collections_tutorials_tutorial_particle_tracing_26_5.png

Part 2: Reflecting boundary conditions#

Next, create a new simulation with 50 time steps and 15 particles in the cylinder.

In addition, enable reflecting boundary conditions in the radial direction. In Struphy this corresponds to logical direction \(\eta_1\):

[15]:

time_opts = Time(dt=0.2, Tend=10.0)
loading_params = LoadingParameters(Np=15, spatial="disc")
boundary_params = BoundaryParameters(bc=("reflect", "periodic", "periodic"))

[16]:

# light-weight model instance
model = Vlasov()

sim = Simulation(model,
                 env=env,
                 time_opts=time_opts,
                 domain=domain,
                 equil=equil,
                 grid=grid,
                 derham_opts=derham_opts,)

[17]:

model.kinetic_ions.set_markers(
    loading_params=loading_params,
    weights_params=weights_params,
    boundary_params=boundary_params,
    saving_params=saving_params,
)

Set propagator options and initial conditions as in the previous section:

[18]:

# propagator options
model.propagators.push_vxb.options = model.propagators.push_vxb.Options()
model.propagators.push_eta.options = model.propagators.push_eta.Options()

# initial conditions (background + perturbation)
perturbation = None
background = maxwellians.Maxwellian3D(n=(1.0, perturbation))

model.kinetic_ions.var.add_background(background)

Run the simulation and plot the resulting trajectories to verify reflections at the radial boundary:

[19]:

out = sim.run()

WARNING: Class "BasisProjectionOperators" called with degree=(1, 1, 1) (interpolation of piece-wise constants should be avoided).

Time stepping: 100%|██████████| 50/50 [00:00<00:00, 359.50step/s]

out.evaluate("<species>/orbits") (or out.orbits.<species>) returns the saved marker orbits as an xarray.Dataset with dimensions t (saved time step) and marker (particle index). Each saved quantity is one variable of this dataset:

  • x, y, z: physical positions,

  • v1, v2, v3: Cartesian velocities (for Particles6D; other particle classes save their own velocity coordinates),

  • weight: time-dependent marker weight.

Each variable carries a description attribute, e.g. orbits.x.attrs["description"]. Use orbits.x.values for a NumPy array of shape (Nt, Np).

[20]:

orbits = out.evaluate("kinetic_ions/orbits")

Nt = orbits.sizes["t"]
Np = orbits.sizes["marker"]


No post-processed data in /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_1, processing with default options (call out.pproc(...) to choose them)

Post-processing path /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_1

No feec fields found in hdf5 file, skipping post-processing of fields.
Evaluation of 15 marker orbits for kinetic_ions

100%|██████████| 51/51 [00:00<00:00, 857.53it/s]

[21]:

import numpy as np

fig = plt.figure()
ax = fig.gca()

colors = ["tab:blue", "tab:orange", "tab:green", "tab:red"]

# create alpha for color scaling
Tend = time_opts.Tend
alpha = np.linspace(1.0, 0.0, Nt + 1)

# loop through particles, plot all time steps
for i in range(Np):
    ax.scatter(orbits.x[:, i], orbits.y[:, i], c=colors[i % 4], alpha=alpha)

circle1 = plt.Circle((0, 0), a2, color="k", fill=False)

ax.add_patch(circle1)
ax.set_aspect("equal")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title(f"{Nt - 1} time steps (full color at t=0)");
../../_images/_collections_tutorials_tutorial_particle_tracing_37_0.png

Part 3: Particles in a cylinder with a magnetic field#

Now extend the setup by adding a static magnetic background field.

You do this by assigning an object from the equils module:

[22]:

B0x = 0.0
B0y = 0.0
B0z = 1.0
equil = equils.HomogenSlab(B0x=B0x, B0y=B0y, B0z=B0z)

To evaluate the equilibrium efficiently in particle kernels, project it onto the spline basis. This requires creating a De Rham complex:

[23]:

bcs = (("free", "free"), None, None)
derham_opts = DerhamOptions(bcs=bcs)

Create a lightweight model instance and configure species options.

Here we choose to remove particles that hit the boundary in radial direction \(\eta_1\):

[24]:

# light-weight model instance
model = Vlasov()

sim_withB = sim.spawn_sister(model=model,
                             equil=equil,
                             derham_opts=derham_opts,)

[25]:

loading_params = LoadingParameters(Np=20)
boundary_params = BoundaryParameters(bc=("remove", "periodic", "periodic"))

model.kinetic_ions.set_markers(
    loading_params=loading_params,
    weights_params=weights_params,
    boundary_params=boundary_params,
    saving_params=saving_params,
)

# propagator options
model.propagators.push_vxb.options = model.propagators.push_vxb.Options()
model.propagators.push_eta.options = model.propagators.push_eta.Options()

# initial conditions (background + perturbation)
perturbation = None
background = maxwellians.Maxwellian3D(n=(1.0, perturbation))

model.kinetic_ions.var.add_background(background)

Run the case, load the data, and plot trajectories to see how the magnetic field modifies particle motion:

[26]:

out_withB = sim_withB.run()

WARNING: Class "BasisProjectionOperators" called with degree=(1, 1, 1) (interpolation of piece-wise constants should be avoided).

Time stepping: 100%|██████████| 50/50 [00:00<00:00, 351.99step/s]

[27]:

orbits = out_withB.evaluate("kinetic_ions/orbits")

Nt = orbits.sizes["t"]
Np = orbits.sizes["marker"]


No post-processed data in /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_1, processing with default options (call out.pproc(...) to choose them)

Post-processing path /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_1

No feec fields found in hdf5 file, skipping post-processing of fields.
Evaluation of 20 marker orbits for kinetic_ions

100%|██████████| 51/51 [00:00<00:00, 828.45it/s]

[28]:

fig = plt.figure()
ax = fig.gca()

colors = ["tab:blue", "tab:orange", "tab:green", "tab:red"]

# create alpha for color scaling
Tend = time_opts.Tend
alpha = np.linspace(1.0, 0.0, Nt + 1)

# loop through particles, plot all time steps
for i in range(Np):
    ax.scatter(orbits.x[:, i], orbits.y[:, i], c=colors[i % 4], alpha=alpha)

circle1 = plt.Circle((0, 0), a2, color="k", fill=False)

ax.add_patch(circle1)
ax.set_aspect("equal")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title(f"{int(Nt - 1)} time steps (full color at t=0)");
../../_images/_collections_tutorials_tutorial_particle_tracing_48_0.png

Part 4: Particles in a tokamak equilibrium#

As a more realistic example, use an ASDEX-Upgrade equilibrium loaded from an EQDSK file.

Instantiate an EQDSKequilibrium with mostly default settings and adjust the density:

[29]:

n1 = 0.0
n2 = 0.0
na = 1.0
equil = equils.EQDSKequilibrium(n1=n1, n2=n2, na=na)

/opt/hostedtoolcache/Python/3.10.21/x64/lib/python3.10/site-packages/struphy/fields_background/equils.py:1694: UserWarning: self.units =<struphy.physics.physics.Units object at 0x7f99748dbac0>, no rescaling performed in EQDSK output.
  warnings.warn(

EQDSKequilibrium inherits from AxisymmMHDequilibrium and CartesianMHDequilibrium, so in principle you could choose different mappings.

For a domain that matches the equilibrium boundary, use the Tokamak mapping:

[30]:

num_elements = (28, 72)
degree = (3, 3)
psi_power = 0.6
psi_shifts = (1e-6, 1.0)
domain = domains.Tokamak(equilibrium=equil, num_elements=num_elements, degree=degree, psi_power=psi_power, psi_shifts=psi_shifts)

The Tokamak domain is a PoloidalSplineTorus. The coordinate relation between Cartesian \((x,y,z)\) and Tokamak \((R,Z,\phi)\) variables is

\[\begin{split}\begin{aligned} x &= R \cos(\phi)\,, \\ y &= -R \sin(\phi)\,, \\ z &= Z\,, \end{aligned}\end{split}\]

where \((R,Z)\) spans the poloidal plane.

Tokamak coordinates are connected to torus coordinates \((r,\theta,\phi)\) through a polar mapping in the poloidal plane:

\[\begin{split}\begin{aligned} R &= R_0 + r \cos(\theta)\,, \\ Z &= r \sin(\theta)\,, \\ \phi &= \phi\,, \end{aligned}\end{split}\]

and torus coordinates map to Struphy logical coordinates \(\boldsymbol\eta=(\eta_1,\eta_2,\eta_3)\in[0,1]^3\) as

\[\begin{split}\begin{aligned} r &= a_1 + (a_2-a_1)\eta_1\,, \\ \theta &= 2\pi\eta_2\,, \\ \phi &= 2\pi\eta_3\,, \end{aligned}\end{split}\]

with \(a_2>a_1\geq 0\) the radial bounds.

You can inspect this relation in mappings such as HollowTorus (more advanced angle parametrizations \(\theta(\eta_1,\eta_2)\) also exist).

As a quick sanity check, plot the magnetic-field magnitude:

  1. in the poloidal plane at \(\phi=0\),

  2. in the top view at \(z=0\).

[31]:

import numpy as np

# logical grid on the unit cube
e1 = np.linspace(0.0, 1.0, 101)
e2 = np.linspace(0.0, 1.0, 101)
e3 = np.linspace(0.0, 1.0, 101)

# move away from the singular point r = 0
e1[0] += 1e-5

[32]:

# logical coordinates of the poloidal plane at phi = 0
eta_poloidal = (e1, e2, 0.0)
# logical coordinates of the top view at theta = 0
eta_topview_1 = (e1, 0.0, e3)
# logical coordinates of the top view at theta = pi
eta_topview_2 = (e1, 0.5, e3)

[33]:

# Cartesian coordinates (squeezed)
x_pol, y_pol, z_pol = domain(*eta_poloidal, squeeze_out=True)
x_top1, y_top1, z_top1 = domain(*eta_topview_1, squeeze_out=True)
x_top2, y_top2, z_top2 = domain(*eta_topview_2, squeeze_out=True)

print(f"{x_pol.shape = }")
print(f"{x_top1.shape = }")
print(f"{x_top2.shape = }")

x_pol.shape = (101, 101)
x_top1.shape = (101, 101)
x_top2.shape = (101, 101)

[34]:

# generate two axes
fig, axs = plt.subplots(2, 1, figsize=(8, 16))
ax = axs[0]
ax_top = axs[1]

# min/max of field strength
equil.domain = domain
Bmax = np.max(equil.absB0(*eta_topview_2, squeeze_out=True))
Bmin = np.min(equil.absB0(*eta_topview_1, squeeze_out=True))
levels = np.linspace(Bmin, Bmax, 51)

# absolute magnetic field at phi = 0
im = ax.contourf(x_pol, z_pol, equil.absB0(*eta_poloidal, squeeze_out=True), levels=levels)

# absolute magnetic field at Z = 0
im_top = ax_top.contourf(x_top1, y_top1, equil.absB0(*eta_topview_1, squeeze_out=True), levels=levels)
ax_top.contourf(x_top2, y_top2, equil.absB0(*eta_topview_2, squeeze_out=True), levels=levels)

# last closed flux surface, poloidal
ax.plot(x_pol[-1], z_pol[-1], color="k")

# last closed flux surface, toroidal
ax_top.plot(x_top1[-1], y_top1[-1], color="k")
ax_top.plot(x_top2[-1], y_top2[-1], color="k")

# limiter, poloidal
ax.plot(equil.limiter_pts_R, equil.limiter_pts_Z, "tab:orange")
ax.axis("equal")
ax.set_xlabel("R")
ax.set_ylabel("Z")
ax.set_title("abs(B) at $\phi=0$")
fig.colorbar(im)
# limiter, toroidal
limiter_Rmax = np.max(equil.limiter_pts_R)
limiter_Rmin = np.min(equil.limiter_pts_R)

thetas = 2 * np.pi * e2
limiter_x_max = limiter_Rmax * np.cos(thetas)
limiter_y_max = -limiter_Rmax * np.sin(thetas)
limiter_x_min = limiter_Rmin * np.cos(thetas)
limiter_y_min = -limiter_Rmin * np.sin(thetas)

ax_top.plot(limiter_x_max, limiter_y_max, "tab:orange")
ax_top.plot(limiter_x_min, limiter_y_min, "tab:orange")
ax_top.axis("equal")
ax_top.set_xlabel("x")
ax_top.set_ylabel("y")
ax_top.set_title("abs(B) at $Z=0$")
fig.colorbar(im_top);
../../_images/_collections_tutorials_tutorial_particle_tracing_57_0.png

As before, build a De Rham complex so the equilibrium can be projected onto the spline basis:

[35]:

num_elements = (32, 72, 1)
grid = grids.TensorProductGrid(num_elements=num_elements)

degree = (3, 3, 1)
bcs = (("free", "free"), None, None)
derham_opts = DerhamOptions(degree=degree, bcs=bcs)

For this example we run 15000 time steps using a second-order splitting scheme:

[36]:

time_opts = Time(dt=0.2, Tend=3000, split_algo="Strang")

Set up a simulation with four hand-picked particle initial conditions to study representative orbit types in this equilibrium:

[37]:

# light-weight model instance
model = Vlasov()

sim_asdex = Simulation(model,
                       env=env,
                       time_opts=time_opts,
                       domain=domain,
                       equil=equil,
                       grid=grid,
                       derham_opts=derham_opts,)

[38]:

# initial particle positions in phase space
initial = (
    (0.501, 0.001, 0.001, 0.0, 0.0450, -0.04),  # co-passing particle
    (0.511, 0.001, 0.001, 0.0, -0.0450, -0.04),  # counter passing particle
    (0.521, 0.001, 0.001, 0.0, 0.0105, -0.04),  # co-trapped particle
    (0.531, 0.001, 0.001, 0.0, -0.0155, -0.04),
)

loading_params = LoadingParameters(Np=4, seed=1608, specific_markers=initial)
boundary_params = BoundaryParameters(bc=("remove", "periodic", "periodic"))

model.kinetic_ions.set_markers(
    loading_params=loading_params,
    weights_params=weights_params,
    boundary_params=boundary_params,
    saving_params=saving_params,
    bufsize=2.0
)

# propagator options
model.propagators.push_vxb.options = model.propagators.push_vxb.Options()
model.propagators.push_eta.options = model.propagators.push_eta.Options()

# initial conditions (background + perturbation)
perturbation = None
background = maxwellians.Maxwellian3D(n=(1.0, perturbation))

model.kinetic_ions.var.add_background(background)

[39]:

out_asdex = sim_asdex.run()

Time stepping: 15001step [00:46, 320.24step/s]

[40]:

orbits = out_asdex.evaluate("kinetic_ions/orbits")

Nt = orbits.sizes["t"]
Np = orbits.sizes["marker"]


No post-processed data in /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_1, processing with default options (call out.pproc(...) to choose them)
/opt/hostedtoolcache/Python/3.10.21/x64/lib/python3.10/site-packages/struphy/fields_background/equils.py:1694: UserWarning: self.units =<struphy.physics.physics.Units object at 0x7f996a07f7f0>, no rescaling performed in EQDSK output.
  warnings.warn(

Post-processing path /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_1

No feec fields found in hdf5 file, skipping post-processing of fields.
Evaluation of 4 marker orbits for kinetic_ions

100%|██████████| 15002/15002 [00:17<00:00, 872.49it/s]

[41]:

import math

colors = ["tab:blue", "tab:orange", "tab:green", "tab:red"]

dt = time_opts.dt
Tend = time_opts.Tend

for i in range(Np):
    r = np.sqrt(orbits.x[:, i] ** 2 + orbits.y[:, i] ** 2)
    # poloidal
    ax.scatter(r, orbits.z[:, i], c=colors[i % 4], s=1)
    # top view
    ax_top.scatter(orbits.x[:, i], orbits.y[:, i], c=colors[i % 4], s=1)

ax.set_title(f"{math.ceil(Tend / dt)} time steps")
ax_top.set_title(f"{math.ceil(Tend / dt)} time steps")
fig

[41]:
../../_images/_collections_tutorials_tutorial_particle_tracing_67_0.png

Part 5: Guiding centers in a tokamak equilibrium#

Finally, repeat the experiment with the guiding-center model and compare the trajectory behavior to full particle tracing:

[42]:

from struphy.models import GuidingCenter

# light-weight model instance
model = GuidingCenter()

[43]:

time_opts = Time(dt=0.1, Tend=100, split_algo="Strang")

sim_gc = sim_asdex.spawn_sister(model=model, time_opts=time_opts)

[44]:

# initial phase space coordinates
# GuidingCenter now uses the Particles5D space (v_para, mu) instead of Particles5Dvperp (v_para, v_perp).
# Convert the perpendicular velocity v_perp used above into the magnetic moment
# mu = v_perp**2 / (2 * |B0|), evaluated at each marker's logical position:
v_perp = 1.72
B0_at_marker = equil.absB0(0.501, 0.001, 0.001, squeeze_out=True)
mu = v_perp**2 / (2 * B0_at_marker)

initial = (
    (0.501, 0.001, 0.001, -1.935, mu),  # co-passing particle
    (0.501, 0.001, 0.001, 1.935, mu),  # counter-passing particle
    (0.501, 0.001, 0.001, -0.6665, mu),  # co-trapped particle
    (0.501, 0.001, 0.001, 0.4515, mu),
)  # counter-trapped particl

loading_params = LoadingParameters(Np=4, seed=1608, specific_markers=initial)
boundary_params = BoundaryParameters(bc=("remove", "periodic", "periodic"))

model.kinetic_ions.set_markers(
    loading_params=loading_params,
    weights_params=weights_params,
    boundary_params=boundary_params,
    saving_params=saving_params,
    bufsize=2.0,
)

# propagator options
model.propagators.push_bxe.options = model.propagators.push_bxe.Options(tol=1e-5)
model.propagators.push_parallel.options = model.propagators.push_parallel.Options(tol=1e-5)

# initial conditions (background + perturbation)
perturbation = None
background = maxwellians.GyroMaxwellian2D(n=(1.0, perturbation), B0=B0_at_marker)

model.kinetic_ions.var.add_background(background)

[45]:

# generate two axes
fig, axs = plt.subplots(2, 1, figsize=(8, 16))
ax = axs[0]
ax_top = axs[1]

# min/max of field strength
equil.domain = domain
Bmax = np.max(equil.absB0(*eta_topview_2, squeeze_out=True))
Bmin = np.min(equil.absB0(*eta_topview_1, squeeze_out=True))
levels = np.linspace(Bmin, Bmax, 51)

# absolute magnetic field at phi = 0
im = ax.contourf(x_pol, z_pol, equil.absB0(*eta_poloidal, squeeze_out=True), levels=levels)

# absolute magnetic field at Z = 0
im_top = ax_top.contourf(x_top1, y_top1, equil.absB0(*eta_topview_1, squeeze_out=True), levels=levels)
ax_top.contourf(x_top2, y_top2, equil.absB0(*eta_topview_2, squeeze_out=True), levels=levels)

# last closed flux surface, poloidal
ax.plot(x_pol[-1], z_pol[-1], color="k")

# last closed flux surface, toroidal
ax_top.plot(x_top1[-1], y_top1[-1], color="k")
ax_top.plot(x_top2[-1], y_top2[-1], color="k")

# limiter, poloidal
ax.plot(equil.limiter_pts_R, equil.limiter_pts_Z, "tab:orange")
ax.axis("equal")
ax.set_xlabel("R")
ax.set_ylabel("Z")
ax.set_title("abs(B) at $\phi=0$")
fig.colorbar(im)
# limiter, toroidal
limiter_Rmax = np.max(equil.limiter_pts_R)
limiter_Rmin = np.min(equil.limiter_pts_R)

thetas = 2 * np.pi * e2
limiter_x_max = limiter_Rmax * np.cos(thetas)
limiter_y_max = -limiter_Rmax * np.sin(thetas)
limiter_x_min = limiter_Rmin * np.cos(thetas)
limiter_y_min = -limiter_Rmin * np.sin(thetas)

ax_top.plot(limiter_x_max, limiter_y_max, "tab:orange")
ax_top.plot(limiter_x_min, limiter_y_min, "tab:orange")
ax_top.axis("equal")
ax_top.set_xlabel("x")
ax_top.set_ylabel("y")
ax_top.set_title("abs(B) at $Z=0$")
fig.colorbar(im_top);
../../_images/_collections_tutorials_tutorial_particle_tracing_72_0.png 
[46]:

out_gc = sim_gc.run()

Time stepping: 1001step [00:04, 222.65step/s]

[47]:

orbits = out_gc.evaluate("kinetic_ions/orbits")

Nt = orbits.sizes["t"]
Np = orbits.sizes["marker"]


No post-processed data in /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_1, processing with default options (call out.pproc(...) to choose them)
/opt/hostedtoolcache/Python/3.10.21/x64/lib/python3.10/site-packages/struphy/fields_background/equils.py:1694: UserWarning: self.units =<struphy.physics.physics.Units object at 0x7f9963e83eb0>, no rescaling performed in EQDSK output.
  warnings.warn(

Post-processing path /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_1

No feec fields found in hdf5 file, skipping post-processing of fields.
Evaluation of 4 marker orbits for kinetic_ions

100%|██████████| 1002/1002 [00:01<00:00, 866.07it/s]

[48]:

import math

colors = ["tab:blue", "tab:orange", "tab:green", "tab:red"]

dt = time_opts.dt
Tend = time_opts.Tend

for i in range(Np):
    r = np.sqrt(orbits.x[:, i] ** 2 + orbits.y[:, i] ** 2)
    # poloidal
    ax.scatter(r, orbits.z[:, i], c=colors[i % 4], s=1)
    # top view
    ax_top.scatter(orbits.x[:, i], orbits.y[:, i], c=colors[i % 4], s=1)

ax.set_title(f"{math.ceil(Tend / dt)} time steps")
ax_top.set_title(f"{math.ceil(Tend / dt)} time steps")
fig

[48]:
../../_images/_collections_tutorials_tutorial_particle_tracing_75_0.png