native parametric tokamak source

docs/source/io_formats/settings.rst

@@ -777,9 +777,9 @@ attributes/sub-elements:
     *Default*: 1.0
 
   :type:
-    Indicator of source type. One of ``independent``, ``file``, ``compiled``, or
-    ``mesh``. The type of the source will be determined by this attribute if it
-    is present.
+    Indicator of source type. One of ``independent``, ``file``, ``compiled``,
+    ``mesh``, or ``tokamak``. The type of the source will be determined by this
+    attribute if it is present.
 
   :particle:
     The source particle type, specified as a PDG number or a string alias (e.g.,
@@ -1015,6 +1015,77 @@ attributes/sub-elements:
     mesh element and follows the format for :ref:`source_element`. The number of
     ``<source>`` sub-elements should correspond to the number of mesh elements.
 
+  For a source with ``type="tokamak"``, the spatial distribution is described by
+  a Miller-style flux-surface parameterization and the following sub-elements
+  are used instead of the ``space`` element:
+
+  :major_radius:
+    The major radius :math:`R_0` of the plasma in [cm].
+
+  :minor_radius:
+    The minor radius :math:`a` of the plasma in [cm]. Must be smaller than
+    ``major_radius``.
+
+  :elongation:
+    The plasma elongation :math:`\kappa` (must be > 0).
+
+  :triangularity:
+    The plasma triangularity :math:`\delta` (must be in [-1, 1]). Negative
+    values describe negative-triangularity plasmas.
+
+  :shafranov_shift:
+    The Shafranov shift :math:`\Delta` in [cm] (must be >= 0 and less than
+    ``minor_radius``/2).
+
+  :r_over_a:
+    A list of normalized minor-radius grid points :math:`r/a`. Must be strictly
+    increasing, start at 0, and end at 1.
+
+  :emission_density:
+    A list of neutron emission densities :math:`S(r)` evaluated at each
+    ``r_over_a`` grid point (arbitrary units, must be non-negative). Only the
+    shape matters, since the profile is normalized internally. Must have the
+    same length as ``r_over_a``.
+
+  :phi_start:
+    The starting toroidal angle in [rad].
+
+    *Default*: 0.0
+
+  :phi_extent:
+    The toroidal angle extent in [rad]. The source is sampled uniformly in
+    :math:`[\phi_\text{start},\ \phi_\text{start} + \phi_\text{extent}]`.
+
+    *Default*: :math:`2\pi`
+
+  :n_alpha:
+    The number of poloidal-angle grid points used to build the sampling CDFs
+    (must be > 2). Larger values reduce discretization bias.
+
+    *Default*: 101
+
+  :vertical_shift:
+    A vertical shift of the plasma center in [cm].
+
+    *Default*: 0.0
+
+  :energy:
+    For a tokamak source, one or more ``energy`` sub-elements specify the
+    neutron energy distribution(s). Either a single distribution is given (used
+    at all radii) or exactly one distribution per ``r_over_a`` grid point is
+    given, in which case the energy is sampled from one of the two
+    distributions bracketing the sampled radius, selected stochastically with
+    probability proportional to the proximity of the radius to each grid point
+    (stochastic interpolation). Each follows the format of a univariate
+    probability distribution (see :ref:`univariate`).
+
+  :time:
+    An optional ``time`` sub-element specifying the time distribution of source
+    particles, following the format of a univariate probability distribution
+    (see :ref:`univariate`).
+
+    *Default*: particles are born at :math:`t=0`
+
   .. note:: Biased sampling can be applied to the spatial and energy distributions
             of a source by using the ``<bias>`` sub-element (see
             :ref:`univariate` for details on how to specify bias distributions).

docs/source/pythonapi/base.rst

@@ -26,6 +26,7 @@ Simulation Settings
    openmc.FileSource
    openmc.CompiledSource
    openmc.MeshSource
+   openmc.TokamakSource
    openmc.SourceParticle
    openmc.VolumeCalculation
    openmc.Settings

docs/source/usersguide/settings.rst

@@ -291,6 +291,51 @@ example, the following would generate a photon source::
 For a full list of all classes related to statistical distributions, see
 :ref:`pythonapi_stats`.
 
+Tokamak Plasma Sources
+----------------------
+
+For fusion applications, the :class:`openmc.TokamakSource` class provides a
+native parametric neutron source for tokamak plasmas. Rather than specifying
+spatial, angular, and energy distributions separately, the source is defined by
+the plasma geometry (using a Miller-style flux-surface parameterization) and a
+radial emission profile. Source sites are sampled directly from the plasma
+volume without rejection.
+
+The plasma shape is described by the major radius :math:`R_0`, minor radius
+:math:`a`, elongation :math:`\kappa`, triangularity :math:`\delta`, and
+Shafranov shift :math:`\Delta`. The neutron birth profile is given as an
+emission density :math:`S(r/a)` tabulated on a normalized minor-radius grid that
+runs from 0 (magnetic axis) to 1 (last closed flux surface); only the shape of
+the profile matters, since it is normalized internally. For example::
+
+  import numpy as np
+
+  r_over_a = np.linspace(0.0, 1.0, 50)
+  emission = (1.0 - r_over_a**2)**2  # peaked on-axis profile
+
+  source = openmc.TokamakSource(
+      major_radius=620.0,    # cm
+      minor_radius=200.0,    # cm
+      elongation=1.8,
+      triangularity=0.45,
+      shafranov_shift=10.0,  # cm
+      r_over_a=r_over_a,
+      emission_density=emission,
+      energy=openmc.stats.muir(e0=14.08e6, m_rat=5.0, kt=20000.0),
+  )
+
+  settings.source = source
+
+The ``energy`` argument accepts either a single
+:class:`~openmc.stats.Univariate` distribution applied at all radii, or a
+sequence with one distribution per ``r_over_a`` grid point to model a
+radially-varying neutron spectrum (energies are then sampled by stochastic
+interpolation between the two distributions bracketing the sampled radius). A
+time distribution can be given with the ``time`` argument; by default,
+particles are born at :math:`t=0`. The toroidal extent can be restricted with
+``phi_start`` and ``phi_extent`` to model a sector of the plasma, and
+``vertical_shift`` translates the plasma center along the z-axis.
+
 File-based Sources
 ------------------
 

include/openmc/math_functions.h

@@ -213,6 +213,20 @@ double exprel(double x);
 //! \return log(1+x)/x without loss of precision near 0
 double log1prel(double x);
 
+//! Evaluate the cylindrical Bessel function of the first kind J_n(x)
+//!
+//! Uses std::cyl_bessel_j where available (e.g., libstdc++). On standard
+//! library implementations lacking the C++17 special math functions (e.g.,
+//! libc++ on Apple Clang/LLVM), falls back to the ascending power series,
+//! which converges to machine precision for the small arguments (|x| <= 2)
+//! used in OpenMC. Unlike std::cyl_bessel_j, negative arguments are handled
+//! via the parity relation J_n(-x) = (-1)^n J_n(x).
+//!
+//! \param n Non-negative integer order of the Bessel function
+//! \param x Real argument
+//! \return J_n(x)
+double cyl_bessel_j(int n, double x);
+
 //! Helper function to get index and interpolation function on an incident
 //! energy grid
 //!

include/openmc/source.h

@@ -7,9 +7,12 @@
 #include <atomic>
 #include <limits>
 #include <unordered_set>
+#include <utility> // for pair
 
 #include "pugixml.hpp"
 
+#include "openmc/array.h"
+#include "openmc/distribution.h"
 #include "openmc/distribution_multi.h"
 #include "openmc/distribution_spatial.h"
 #include "openmc/memory.h"
@@ -260,6 +263,140 @@ class MeshSource : public Source {
   vector<unique_ptr<IndependentSource>> sources_; //!< Source distributions
 };
 
+//==============================================================================
+//! Parametric tokamak plasma neutron source
+//!
+//! This source samples neutron positions from a tokamak plasma geometry using
+//! Miller-style flux surface parameterization with user-specified emission
+//! profiles and energy distributions.
+//!
+//! Flux surface parameterization:
+//!   R = R0 + r*cos(alpha + delta*sin(alpha)) + Delta*(1 - (r/a)^2)
+//!   Z = kappa * r * sin(alpha)
+//!
+//! The sampling algorithm:
+//! 1. Sample minor radius r from precomputed CDF of S(r) * Jacobian
+//! 2. Sample poloidal angle alpha from conditional P(alpha|r) using mixture
+//!    of precomputed CDFs weighted by functions of r
+//! 3. Sample energy and time from user-provided distribution(s)
+//! 4. Sample isotropic direction
+//! 5. Sample toroidal angle phi uniformly in [phi_start, phi_start +
+//! phi_extent]
+//! 6. Transform (r, alpha, phi) to Cartesian (x, y, z), applying the optional
+//!    vertical shift of the plasma center
+//!
+//! The user provides the emission density S(r) directly (e.g., from transport
+//! codes like TRANSP, ASTRA, etc.), allowing full flexibility in reaction
+//! physics calculations. S(r) is a profile in arbitrary units sampled on the
+//! r_over_a grid; only its shape matters, since it is normalized internally.
+//! Energy distributions can be specified as either a single distribution for
+//! all r, or one distribution per radial point.
+//==============================================================================
+
+class TokamakSource : public Source {
+public:
+  // Constructors
+  explicit TokamakSource(pugi::xml_node node);
+
+  //! Sample from the tokamak source distribution
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Sampled site
+  SourceSite sample(uint64_t* seed) const override;
+
+private:
+  //==========================================================================
+  // Private methods
+
+  //! Precompute data structures for efficient sampling
+  void precompute_sampling_cdfs();
+
+  //! Sample minor radius from marginal CDF
+  //! \param seed Pseudorandom seed pointer
+  //! \return Sampled r/a value
+  double sample_r_over_a(uint64_t* seed) const;
+
+  //! Sample poloidal angle given r using mixture of precomputed CDFs
+  //! \param r_norm Normalized minor radius r/a
+  //! \param seed Pseudorandom seed pointer
+  //! \return Sampled poloidal angle alpha [rad]
+  double sample_poloidal_angle(double r_norm, uint64_t* seed) const;
+
+  //! Compute the k-th mixture weight w_k(r) * I_hat_k for poloidal sampling
+  //! \param k Basis function index (0-5)
+  //! \param r Normalized minor radius r/a
+  //! \return Mixture weight for component k
+  double mixture_weight(int k, double r) const;
+
+  //! Sample energy from the distribution(s)
+  //! \param r_norm Normalized minor radius r/a (for distribution selection)
+  //! \param seed Pseudorandom seed pointer
+  //! \return (Sampled energy [eV], importance weight)
+  std::pair<double, double> sample_energy(double r_norm, uint64_t* seed) const;
+
+  //! Transform from flux coordinates (r, alpha, phi) to Cartesian (x, y, z)
+  //! \param r Minor radius [cm]
+  //! \param alpha Poloidal angle [rad]
+  //! \param phi Toroidal angle [rad]
+  //! \return Position in Cartesian coordinates [cm]
+  Position flux_to_cartesian(double r, double alpha, double phi) const;
+
+  //==========================================================================
+  // Data members
+
+  // Emission profile (input)
+  vector<double> r_over_a_;         //!< Normalized minor radius grid points
+  vector<double> emission_density_; //!< Emission density S(r) at grid points
+
+  // Energy distribution(s): either 1 for all r, or one per r point
+  vector<unique_ptr<Distribution>> energy_dists_;
+
+  // Time distribution (defaults to a delta distribution at t=0)
+  UPtrDist time_;
+
+  // Angular distribution (isotropic)
+  UPtrAngle angle_;
+
+  // Tokamak geometry parameters
+  double major_radius_;    //!< Major radius R0 [cm]
+  double minor_radius_;    //!< Minor radius a [cm]
+  double elongation_;      //!< Elongation kappa
+  double triangularity_;   //!< Triangularity delta
+  double shafranov_shift_; //!< Shafranov shift Delta [cm]
+  double vertical_shift_;  //!< Vertical shift of plasma center [cm]
+
+  // Normalized geometry parameters (precomputed for efficiency)
+  double epsilon_;     //!< Inverse aspect ratio a/R0
+  double delta_tilde_; //!< Normalized Shafranov shift Delta/a
+
+  // Toroidal angle bounds
+  double phi_start_;  //!< Starting toroidal angle [rad]
+  double phi_extent_; //!< Toroidal angle extent [rad]
+
+  // Precomputed data for radial sampling
+  vector<double> radial_cdf_; //!< CDF for radial sampling (on r_over_a_ grid)
+
+  // Coefficients of the radial geometric polynomial: A*r - B*r^2 - C*r^3
+  // Also used as the analytical normalization for poloidal mixture weights
+  double radial_poly_a_; //!< 1 + ε·Δ̃
+  double radial_poly_b_; //!< (3/8)·c₁·ε
+  double radial_poly_c_; //!< 2·ε·Δ̃
+
+  // Precomputed Bernstein basis functions for poloidal angle sampling.
+  // Using the factorization f(r_tilde, alpha) = R_tilde x J_tilde where:
+  //   R_tilde = b0*(1-r)^2 + 2*b1*r*(1-r) + b2*r^2  (Bernstein quadratic)
+  //   J_tilde = b3*(1-r) + b4*r                      (Bernstein linear)
+  // The product gives 6 non-negative basis functions g_k(alpha) with
+  // weights w_k(r_tilde) that are products of Bernstein polynomials.
+  // CDFs are computed on [0, pi] exploiting up-down symmetry.
+  static constexpr int N_POLOIDAL_BASIS = 6; //!< Number of basis functions
+  int n_alpha_;                        //!< Number of poloidal angle grid points
+  vector<double> poloidal_alpha_grid_; //!< alpha grid for CDF lookup [0, pi]
+  array<vector<double>, N_POLOIDAL_BASIS>
+    poloidal_cdfs_; //!< CDFs for each basis function g_k(alpha)
+  array<double, N_POLOIDAL_BASIS>
+    poloidal_integrals_; //!< Integrals of g_k(alpha) over [0, pi]
+};
+
 //==============================================================================
 // Functions
 //==============================================================================