Tracing markers backwards in time

This example shows how to trace markers backwards in time.

For starters it is worth emphasizing that this tutorial is only for tracing marker orbits backwards in time and not for reverting the physical process. Activating the time reversal option, REVERSE_TIME=1, makes the time run backwards in simulations so that one can e.g. see the orbit at which a marker hit the wall or FILD. However, it does not reverse collisional or other stochastic process so these must be disabled. Use BMC for proper backwards Monte Carlo modelling that include these processess and can properly run a slowing-down simulation in reverse.

First set up the test data.

import numpy as np
import matplotlib.pyplot as plt
from a5py import Ascot

a5 = Ascot("ascot.h5", create=True)
a5.data.create_input("bfield analytical iter circular")
a5.data.create_input("plasma flat")
a5.data.create_input("E_TC")
a5.data.create_input("N0_1D")
a5.data.create_input("Boozer")
a5.data.create_input("MHD_STAT")
a5.data.create_input("asigma_loc")

print("Inputs created")

Inputs created

As for the wall, any wall data would do but we wish to do some nice visualizations in 3D so we choose a 3D wall for this exercise.

from a5py.ascot5io.wall import wall_3D

rad  = 2.0
pol  = np.linspace(0, 2*np.pi, 181)[:-1]
w2d = {"nelements":180,
       "r":6.4 + rad*np.cos(pol), "z":rad*np.sin(pol)}
w3d = wall_3D.convert_wall_2D(180, **w2d)
a5.data.create_input("wall_3D", **w3d, desc="REVOLVED")

'wall_3D_3502812368'

Now let’s create a few markers which we trace for a few moments as if we were modelling the slowing-down process. This is an ordinary forward in time simulation, we will use the results to run the reversed simulation (although there is nothing preventing from running the reversed simulation without the forward simulation).

from a5py.ascot5io.marker import Marker
nmrk = 10
mrk = Marker.generate("gc", n=nmrk, species="alpha")
mrk["energy"][:] = 3.5e6
mrk["pitch"][:]  = 0.6
mrk["r"][:]      = np.linspace(7.6,8.0,nmrk)
a5.data.create_input("gc", **mrk, activate=True)

from a5py.ascot5io.options import Opt

opt = Opt.get_default()
opt.update({
    # Reversed time does not work well with the adaptive step at the moment so disable it
    "SIM_MODE":2, "ENABLE_ADAPTIVE":0,
    "ENDCOND_SIMTIMELIM":1, "ENDCOND_LIM_SIMTIME":1.0e-4, "ENDCOND_WALLHIT":1,
    "ENABLE_ORBIT_FOLLOWING":1, "ENABLE_COULOMB_COLLISIONS":1,
    "ENABLE_ORBITWRITE":1, "ORBITWRITE_MODE":1, "ORBITWRITE_INTERVAL":0,
    "ORBITWRITE_NPOINT":10**4,
})
a5.data.create_input("opt", **opt, desc="TIMEFORWARD", activate=True)

'opt_0365348477'

Run the simulation

import subprocess
subprocess.run(["./../../build/ascot5_main", "--d=\"GREATESTHITSVOL2\""])
print("Simulation completed")

ASCOT5_MAIN
Tag e389e6c
Branch docs

Initialized MPI, rank 0, size 1.

Reading and initializing input.

Input file is ascot.h5.

Reading options input.
Active QID is 0365348477
Options read and initialized.

Reading magnetic field input.
Active QID is 1860762275

Analytical tokamak magnetic field (B_GS)
Psi at magnetic axis (6.618 m, -0.000 m)
-5.410 (evaluated)
-5.410 (given)
Magnetic field on axis:
B_R = 0.000 B_phi = 4.965 B_z = -0.000
Number of toroidal field coils 0
Estimated memory usage 0.0 MB
Magnetic field read and initialized.

Reading electric field input.
Active QID is 0595877904

Trivial Cartesian electric field (E_TC)
E_x = 0.000000e+00, E_y = 0.000000e+00, E_z = 0.000000e+00
Estimated memory usage 0.0 MB
Electric field read and initialized.

Reading plasma input.
Active QID is 0546981494

1D plasma profiles (P_1D)
Min rho = 0.00e+00, Max rho = 1.00e+01, Number of rho grid points = 100, Number of ion species = 1
Species Z/A  charge [e]/mass [amu] Density [m^-3] at Min/Max rho    Temperature [eV] at Min/Max rho
   1  /  1     1  /  1.000             1.00e+21/1.00e+00                1.00e+04/1.00e+04
[electrons]   -1  /  0.001             1.00e+21/1.00e+00                1.00e+04/1.00e+04
Quasi-neutrality is (electron / ion charge density) 1.00
Estimated memory usage 0.0 MB
Plasma data read and initialized.

Reading neutral input.
Active QID is 4138853207

1D neutral density and temperature (N0_1D)
Grid:  nrho =    3   rhomin = 0.000   rhomax = 2.000
 Number of neutral species = 1
Species Z/A   (Maxwellian)
        1/  1 (1)
Estimated memory usage 0.0 MB
Neutral data read and initialized.

Reading wall input.
Active QID is 3502812368

3D wall model (wall_3D)
Number of wall elements 64800
Spanning xmin = -8.500 m, xmax = 8.500 m
         ymin = -8.500 m, ymax = 8.500 m
         zmin = -2.100 m, zmax = 2.100 m
Estimated memory usage 4.4 MB
Wall data read and initialized.

Reading boozer input.
Active QID is 4217749134

Boozer input
psi grid: n =    6 min = 0.000 max = 1.000
thetageo grid: n =   18
thetabzr grid: n =   10
Boozer data read and initialized.

Reading MHD input.
Active QID is 3055972933

MHD (stationary) input
Grid: nrho =    6 rhomin = 0.000 rhomax = 1.000

Modes:
(n,m) = ( 1, 3) Amplitude = 0.1 Frequency =   1 Phase =   0
(n,m) = ( 2, 4) Amplitude =   2 Frequency = 1.5 Phase = 0.785
Estimated memory usage 0.0 MB
MHD data read and initialized.

Reading atomic reaction input.
Active QID is 1110749216

Found data for 1 atomic reactions:
Reaction number / Total number of reactions = 1 / 1
  Reactant species Z_1 / A_1, Z_2 / A_2     = 1 / 1, 1 / 1
  Min/Max energy                            = 1.00e+03 / 1.00e+04
  Min/Max density                           = 1.00e+18 / 1.00e+20
  Min/Max temperature                       = 1.00e+03 / 1.00e+04
  Number of energy grid points              = 3
  Number of density grid points             = 4
  Number of temperature grid points         = 5
Estimated memory usage 0.0 MB
Atomic reaction data read and initialized.

Reading marker input.
Active QID is 2692907280

Loaded 10 guiding centers.
Marker data read and initialized.

All input read and initialized.

Initializing marker states.
Estimated memory usage 0.0 MB.
Marker states initialized.
Initialized diagnostics, 12.2 MB.

Preparing output.
Note: Output file ascot.h5 is already present.

The qid of this run is 0113848382

Inistate written.
Simulation begins; 4 threads.
Simulation complete.
Simulation finished in 0.775830 s
Endstate written.

Combining and writing diagnostics.

Writing diagnostics output.
Writing orbit diagnostics.

Diagnostics output written.
Diagnostics written.

Summary of results:
        7 markers had end condition Sim time limit
        3 markers had end condition Wall collision

          No markers were aborted.

Done.
Simulation completed

Few markers hit the wall. We select one of them and run the orbit backwards in time.

a5 = Ascot("ascot.h5")

# Pick one of the markers that hit the wall
ids = a5.data.active.getstate("ids", endcond="wall")[0]

# When converting endstate to marker input, the input type and the simulation
# mode that was used should preferably match
mrk = a5.data.active.getstate_markers("gc", ids=ids)

opt = a5.data.active.options.read()
opt.update({
    # Reverse time and turn off the collisions
    "REVERSE_TIME":1, "ENABLE_COULOMB_COLLISIONS":0,
    # When using the reversed time, the simulation is stopped if time < LIM_SIMTIME.
    # To follow marker back to its original position, we set LIM_SIMTIME=0 which was
    # the marker time at the start of the forward simulation.
    # Note that MAX_MILEAGE works same way both in forward and reversed simulations.
    "ENDCOND_LIM_SIMTIME":0.0,
    "ENDCOND_WALLHIT":0
})

Since we are only tracing a single marker, and we are only interested in its orbit, the reversed simulation can be done conveniently using the virtual run (but ordinary run would work as well).

Run the simulation and plot both the forward orbit and the reversed orbit.

a5.simulation_initinputs()
a5.simulation_initmarkers(**mrk)
a5.simulation_initoptions(**opt)
vrun = a5.simulation_run(printsummary=False)

fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)

a5.data.active.plotorbit_trajectory("r", "z", ids=mrk["ids"], axes=ax1)
a5.input_plotwallcontour(axes=ax1)

vrun.plotorbit_trajectory("r", "z", axes=ax2)
a5.input_plotwallcontour(axes=ax2)

plt.show()

This is how time reversed simulations are done. They can be used to map the phase-space region where losses originate (for collisionless losses) or to make nice visualizations of marker orbits when they hit the wall or FILD.

fig, ax = plt.subplots()
vrun.plotwall_3dstill(cpos=(-7.0,0.0,0.5), cfoc=(-3.8,-3.5,1.5), cang=(0,-20,-90), orbit=mrk["ids"], axes=ax)
a5.simulation_free()

/usr/share/miniconda/envs/ascot-dev/lib/python3.10/site-packages/pyvista/core/pointset.py:1365: PyVistaDeprecationWarning: The current behavior of `pv.PolyData.n_faces` has been deprecated.
                Use `pv.PolyData.n_cells` or `pv.PolyData.n_faces_strict` instead.
                See the documentation in '`pv.PolyData.n_faces` for more information.
  warnings.warn(
/usr/share/miniconda/envs/ascot-dev/lib/python3.10/site-packages/pyvista/jupyter/notebook.py:37: UserWarning: Failed to use notebook backend:

No module named 'trame'

Falling back to a static output.
  warnings.warn(