Creating fusion source with AFSI
This example shows how to create a fusion source with AFSI in various cases.
AFSI (Ascot Fusion Source Integrator) calculates the 5D distributions of fusion reaction products based on the 5D distributions of the reactants. The algorithm iterates over all cells in cylindrical coordinates, and on each location creates fusion products via Monte Carlo sampling. First, a pair of markers is sampled from the reactant (velocity) distributions and the velocities of the fusion products for a given reaction are calculated using the initial velocities. The reaction probability, which depends on the relative velocity between the reactants, is used to weight the products which in turn are used to calculate the 5D distributions for the output. Several pairs are created before the algorithm proceeds to the next cylindrical cell.
The Monte Carlo sampling makes AFSI capable of not only calculating reactions for thermal (Maxwellian) species but for arbitrary distributions. In practice this often means NBI slowing-down distributions and, hence, the arbitrary distribution is referred to as beam when working with AFSI.
Note that AFSI calculates the distributions for both products, i.e., one can use AFSI to calculate the neutron source as well.
AFSI has three modes of operation based on what the inputs are:
thermalcalculates the fusion products for two Maxwellian species. The population can also react itself as is the case with D-D fusion.beam-thermalcalculates the fusion products for a Maxwellian species that is interacting with an arbitrary distribution.beam-beamcalculates the fusion products for two arbitrary distributions. A distribution can interact with itself, e.g. if there is a single deuterium beam.
AFSI requires only the magnetic field and plasma data for it to work, but we also need to get the beam ion distribution somewhere for this tutorial. The following cell initializes the inputs and runs a short slowing-down simulation for a deuterium beam. Its contents are not relevant for this tutorial.
[1]:
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", desc="ANALYTICAL")
# DT-plasma
nrho = 101
rho = np.linspace(0, 2, nrho).T
prof = (rho<=1)*(1.0 - rho**(3.0/2))**3 + 1e-6
edens = 2e21 * prof
etemp = 1e4 * np.ones((nrho, 1))
idens = 1e21 * np.tile(prof,(2,1)).T
itemp = 1e4 * np.ones((nrho, 1))
edens[rho>=1] = 1
idens[rho>=1,:] = 1
pls = {
"nrho" : nrho, "nion" : 2, "rho" : rho,
"anum" : np.array([2, 3]), "znum" : np.array([1, 1]),
"mass" : np.array([2.014, 3.016]), "charge" : np.array([1, 1]),
"edensity" : edens, "etemperature" : etemp,
"idensity" : idens, "itemperature" : itemp}
a5.data.create_input("plasma_1D", **pls, desc="DT")
# These inputs are not relevant for AFSI
a5.data.create_input("wall rectangular")
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")
# Generate deuterium markers with Ekin=1MeV and run a short slowing down simulation.
from a5py.ascot5io.marker import Marker
mrk = Marker.generate("gc", n=100, species="deuterium")
mrk["energy"][:] = 1.0e6
mrk["pitch"][:] = 0.99 - 0.2 * np.random.rand(100,)
mrk["r"][:] = np.linspace(6.2, 7.2, 100)
a5.data.create_input("gc", **mrk)
from a5py.ascot5io.options import Opt
opt = Opt.get_default()
opt.update({
"SIM_MODE":2, "ENABLE_ADAPTIVE":1,
"ENDCOND_SIMTIMELIM":1, "ENDCOND_MAX_MILEAGE":0.05,
"ENABLE_ORBIT_FOLLOWING":1, "ENABLE_COULOMB_COLLISIONS":1,
"ENABLE_DIST_5D":1,
"DIST_MIN_R":4.3, "DIST_MAX_R":8.3, "DIST_NBIN_R":50,
"DIST_MIN_PHI":0, "DIST_MAX_PHI":360, "DIST_NBIN_PHI":1,
"DIST_MIN_Z":-2.0, "DIST_MAX_Z":2.0, "DIST_NBIN_Z":50,
"DIST_MIN_PPA":-1.3e-19, "DIST_MAX_PPA":1.3e-19, "DIST_NBIN_PPA":100,
"DIST_MIN_PPE":0, "DIST_MAX_PPE":1.3e-19, "DIST_NBIN_PPE":50,
"DIST_MIN_TIME":0, "DIST_MAX_TIME":1.0, "DIST_NBIN_TIME":1,
"DIST_MIN_CHARGE":1, "DIST_MAX_CHARGE":3, "DIST_NBIN_CHARGE":1,
})
a5.data.create_input("opt", **opt, desc="SLOWINGDOWN")
import subprocess
subprocess.run(["./../../build/ascot5_main", "--d=\"BEAMSD\""])
print("Simulation complete")
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 2886371688
Options read and initialized.
Reading magnetic field input.
Active QID is 0999871598
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 1388046862
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 0009326458
1D plasma profiles (P_1D)
Min rho = 0.00e+00, Max rho = 2.00e+00, Number of rho grid points = 101, Number of ion species = 2
Species Z/A charge [e]/mass [amu] Density [m^-3] at Min/Max rho Temperature [eV] at Min/Max rho
1 / 2 1 / 2.014 1.00e+21/1.00e+00 1.00e+04/1.00e+04
1 / 3 1 / 3.016 1.00e+21/1.00e+00 1.00e+04/1.00e+04
[electrons] -1 / 0.001 2.00e+21/1.00e+00 1.00e+04/1.00e+04
Quasi-neutrality is (electron / ion charge density) 2.00
Estimated memory usage 0.0 MB
Plasma data read and initialized.
Reading neutral input.
Active QID is 3647474105
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 0881241138
2D wall model (wall_2D)
Number of wall elements = 20, R extend = [4.10, 8.40], z extend = [-3.90, 3.90]
Estimated memory usage 0.0 MB
Wall data read and initialized.
Reading boozer input.
Active QID is 2966013760
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 3651373014
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 2709402688
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 3211941356
Loaded 100 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, 95.4 MB.
Preparing output.
Note: Output file ascot.h5 is already present.
The qid of this run is 1525920612
Inistate written.
Simulation begins; 4 threads.
Simulation complete.
Simulation finished in 10.976247 s
Endstate written.
Combining and writing diagnostics.
Writing diagnostics output.
Writing 5D distribution.
Diagnostics output written.
Diagnostics written.
Summary of results:
100 markers had end condition Sim time limit
No markers were aborted.
Done.
Simulation complete
AFSI is used via the Ascot object using the attribute afsi which has all the relevant methods. To use AFSI, we need to specify what reaction is that we are considering. To list possible reactions:
[2]:
a5.afsi.reactions()
[2]:
dict_keys(['DT_He4n', 'DHe3_He4p', 'DD_Tp', 'DD_He3n'])
To run AFSI on thermal mode, we need to specify \((R,z)\), and optionally \(\phi\), grid on which the distributions are calculated. Since both reactant populations are Maxwellian, there is no need to define velocity grids for those as AFSI can use the analytical distribution. Specifying the velocity grid for the products is required however.
Once inputs are set, we can call thermal function that computes the product distributions for us.
[3]:
# Grid spanning the whole plasma
rmin = 4.2; rmax = 8.2; nr = 50
zmin = -2.0; zmax = 2.0; nz = 50
# Run AFSI which creates a fusion alpha distribution
a5.afsi.thermal(
"DT_He4n",
rmin, rmax, nr, zmin, zmax, nz,
minphi=0, maxphi=360, nphi=1, nmc=1000,
minppara=-1.3e-19, maxppara=1.3e-19, nppara=80,
minpperp=0, maxpperp=1.3e-19, npperp=40)
[3]:
(<a5py.ascotpy.ascot2py.struct_c__SA_dist_5D_offload_data at 0x7f47adaeecc0>,
<a5py.ascotpy.ascot2py.struct_c__SA_dist_5D_offload_data at 0x7f47adaed8c0>)
The result is stored in the HDF5 file in a similar run group (found in a5.data) as the ASCOT5 simulations are, and the group works in a similar manner. The group contains information on the magnetic field and plasma inputs that were used and provides methods to access the product distributions. However, now the distributions are called prod1 and prod2 instead of 5d, but they act exactly as a 5d distribution would.
[4]:
# Convert to E-xi as then we can make dist more compact in momentum space
alphadist = a5.data.active.getdist("prod1", exi=True, ekin_edges=20, pitch_edges=10)
alphadist.integrate(time=np.s_[:], charge=np.s_[:]) # Distribution must be 5D
rzdist = alphadist.integrate(copy=True, phi=np.s_[:], ekin=np.s_[:], pitch=np.s_[:])
exidist = alphadist.integrate(copy=True, phi=np.s_[:], r=np.s_[:], z=np.s_[:])
fig = plt.figure(figsize=(10,5))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
rzdist.plot(axes=ax1)
exidist.plot(axes=ax2)
plt.show()
Next, we run the beam-thermal fusion between the deuterium beam (whose distribution was calculated earlier) and thermal deuterium. This time we use beamthermal function and provide the beam 5D function as an input. Note that this time we can’t specify a cylindrical grid for the outputs, as the same grid is used as in the beam distribution.
The results is again stored as a new run from which the data is read as usual. Take a note how the product distribution is strongly anisotropic; the product velocities consists of both the kinetic energy released in the reaction and the (beam) velocity that the reactants had.
[5]:
beamdist = a5.data.BEAMSD.getdist("5d")
a5.afsi.beamthermal("DD_Tp", beamdist, swap=True)
alphadist = a5.data.active.getdist("prod1", exi=True, ekin_edges=20, pitch_edges=10)
rzdist = alphadist.integrate(copy=True, phi=np.s_[:], ekin=np.s_[:], pitch=np.s_[:])
exidist = alphadist.integrate(copy=True, phi=np.s_[:], r=np.s_[:], z=np.s_[:])
fig = plt.figure(figsize=(10,5))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
rzdist.plot(axes=ax1)
exidist.plot(axes=ax2)
plt.show()
AFSI5
Tag e389e6c
Branch docs
Note: Output file /home/runner/work/ascot5/ascot5/doc/tutorials/ascot.h5 is already present.
The qid of this run is 1270549620
Done
Finally, we demonstrate the beam-beam fusion. Providing the beambeam function with just a single distribution calculates the fusion product for the beam interacting with itself. The other beam would be provided with beam2 argument (which must have the same abscissae as the first beam distribution).
[6]:
a5 = Ascot("ascot.h5")
beamdist = a5.data.BEAMSD.getdist("5d")
a5.afsi.beambeam("DD_Tp", beamdist)
alphadist = a5.data.active.getdist("prod1", exi=True, ekin_edges=20, pitch_edges=10)
alphadist.integrate(time=np.s_[:], charge=np.s_[:]) # Distribution must be 5D
rzdist = alphadist.integrate(copy=True, phi=np.s_[:], ekin=np.s_[:], pitch=np.s_[:])
exidist = alphadist.integrate(copy=True, phi=np.s_[:], r=np.s_[:], z=np.s_[:])
fig = plt.figure(figsize=(10,5))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
rzdist.plot(axes=ax1)
exidist.plot(axes=ax2)
plt.show()
AFSI5
Tag e389e6c
Branch docs
Note: Output file /home/runner/work/ascot5/ascot5/doc/tutorials/ascot.h5 is already present.
The qid of this run is 0489342356
Done
