From f4d471690007440dd4efb8951eccf9e7f77d1420 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 10 Oct 2026 14:09:14 -0500 Subject: [PATCH 01/15] Sample stellarator angles with radius-dependent cell bounds --- docs/source/usersguide/settings.rst | 14 ++ include/openmc/source.h | 7 +- src/source.cpp | 190 ++++++++++++++++---- tests/unit_tests/test_source_stellarator.py | 42 +++++ 4 files changed, 219 insertions(+), 34 deletions(-) diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index c793f815955..bf2ab4ebbb6 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -476,6 +476,20 @@ radial proposals are resampled within the same interval. The two angles are sampled from the conditional distribution :math:`p(\theta,\zeta|\rho) \propto R\,|\tau|` (with :math:`\tau` the poloidal-plane Jacobian) by rejection. +Angular proposals use a 32 by 32 partition of one field period. Each cell has +a bound that varies quadratically across the radial interval. OpenMC selects +cells in proportion to these bounds, samples uniformly inside the chosen cell, +and accepts with probability equal to the density divided by the cell bound. +Rejected angular proposals retain the radius and select a new cell. A uniformly +chosen field period maps the accepted point to the full torus. + +The bounds are constructed during C++ initialization from the density's Fourier +coefficients, using a Taylor bound on each cell. This requires no additional +dependencies or Python preprocessing. Sampling checks for a density exceeding +its bound and reports an error if one is encountered. Bounds and alias tables +require approximately 72 KiB per radial interval on a 64-bit platform, in +addition to the geometry coefficients and temporary initialization arrays. + Sources are most conveniently created directly from equilibrium code output:: # From a VMEC wout NetCDF file diff --git a/include/openmc/source.h b/include/openmc/source.h index 11079a07067..c35a186b264 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -17,6 +17,7 @@ #include "openmc/distribution_spatial.h" #include "openmc/memory.h" #include "openmc/particle_type.h" +#include "openmc/tensor.h" #include "openmc/vector.h" namespace openmc { @@ -509,8 +510,10 @@ class StellaratorSource : public Source { vector> radial_pdf_; DiscreteIndex radial_bins_; - // Per-radial-bin majorant of R*tau for rejection sampling (n_rho - 1) - vector envelope_; + // Angular cell bounds for the three quadratic Bernstein components per bin. + tensor::Tensor angular_bounds_; + vector> angular_cells_; + vector> angular_totals_; // Global sign of the flux-coordinate Jacobian (+1 or -1) double jacobian_sign_ {1.0}; diff --git a/src/source.cpp b/src/source.cpp index 0fee2678717..55b79c4c438 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -6,8 +6,9 @@ #include // for max #include // for sin, cos, abs -#include // for accumulate -#include // for move +#include +#include // for accumulate +#include // for move #ifdef HAS_DYNAMIC_LINKING #include // for dlopen, dlsym, dlclose, dlerror @@ -1176,23 +1177,105 @@ SourceSite TokamakSource::sample(uint64_t* seed) const namespace { -//! Maximum over t in [0, 1] of the quadratic through q(0)=f0, q(1/2)=fh, -//! q(1)=f1. Exact, since the sampled density is quadratic in the radial -//! interpolation parameter within a bin. -double quad_max01(double f0, double fh, double f1) +constexpr int ANGULAR_BINS = 32; +using Complex = std::complex; + +// Separable Fourier sums. The small transforms used during source setup do +// not require an FFT dependency. Each matrix is indexed [output, input]. +template +vector angular_transform(const vector& values, int ni, int nj, + const vector& left, const vector& right) { - double c = 2.0 * (f0 + f1 - 2.0 * fh); - double b = f1 - f0 - c; - double m = std::max(f0, f1); - if (c < 0.0) { - double t_star = -b / (2.0 * c); - if (t_star > 0.0 && t_star < 1.0) { - m = std::max(m, f0 + t_star * (b + t_star * c)); + int no = left.size() / ni; + int np = right.size() / nj; + vector intermediate(no * nj), result(no * np); + for (int o = 0; o < no; ++o) { + for (int i = 0; i < ni; ++i) { + auto a = left[o * ni + i]; + for (int j = 0; j < nj; ++j) + intermediate[o * nj + j] += a * values[i * nj + j]; + } + for (int p = 0; p < np; ++p) { + for (int j = 0; j < nj; ++j) + result[o * np + p] += intermediate[o * nj + j] * right[p * nj + j]; } } - return m; + return result; } +// Bound a real Fourier polynomial on 32x32 angular cells. Recover the density +// coefficients from its resolved quadrature grid, retaining cancellations in +// the Jacobian. Taylor's theorem bounds each cell using its center value, +// gradient, and a Fourier bound on the second directional derivative. +class AngularBounds { +public: + AngularBounds(int nt, int nz, int m, int n) + : nt_ {nt}, nz_ {nz}, m_ {m}, n_ {n} + { + auto forward = [](int points, int degree) { + vector matrix((2 * degree + 1) * points); + for (int k = -degree; k <= degree; ++k) + for (int j = 0; j < points; ++j) + matrix[(k + degree) * points + j] = + std::polar(1.0 / points, -2 * PI * k * j / points); + return matrix; + }; + theta_ = forward(nt, m); + zeta_ = forward(nz, n); + inverse_.resize(ANGULAR_BINS * ANGULAR_BINS); + for (int i = 0; i < ANGULAR_BINS; ++i) + for (int j = 0; j < ANGULAR_BINS; ++j) + inverse_[i * ANGULAR_BINS + j] = + std::polar(1.0, 2 * PI * i * j / ANGULAR_BINS); + } + + void compute( + const vector& values, double* upper, double* negative_upper) const + { + auto coeff = angular_transform(values, nt_, nz_, theta_, zeta_); + array, 3> folded; + for (auto& v : folded) + v.resize(ANGULAR_BINS * ANGULAR_BINS); + double norm = 0.0, curvature = 0.0; + double h = PI / ANGULAR_BINS; + for (int m = -m_; m <= m_; ++m) { + for (int n = -n_; n <= n_; ++n) { + auto c = coeff[(m + m_) * (2 * n_ + 1) + n + n_]; + double k = std::abs(m) + std::abs(n); + norm += std::abs(c); + curvature += std::abs(c) * k * k; + // Fold frequencies onto the cell-center grid, preserving their phase + // and derivative factors even when they exceed its Nyquist frequency. + int i = (m % ANGULAR_BINS + ANGULAR_BINS) % ANGULAR_BINS; + int j = (n % ANGULAR_BINS + ANGULAR_BINS) % ANGULAR_BINS; + c *= std::polar(1.0, (m + n) * h); + folded[0][i * ANGULAR_BINS + j] += c; + folded[1][i * ANGULAR_BINS + j] += Complex(0, m) * c; + folded[2][i * ANGULAR_BINS + j] += Complex(0, n) * c; + } + } + for (auto& v : folded) + v = angular_transform(v, ANGULAR_BINS, ANGULAR_BINS, inverse_, inverse_); + // Small allowance for rounding in the Fourier sums and density evaluation. + double remainder = 0.5 * h * h * curvature + 1.e-10 * norm; + for (int c = 0; c < ANGULAR_BINS * ANGULAR_BINS; ++c) { + double excursion = + h * (std::abs(folded[1][c].real()) + std::abs(folded[2][c].real())) + + remainder; + double center = folded[0][c].real(); + if (!std::isfinite(center) || !std::isfinite(excursion)) + fatal_error("StellaratorSource: nonfinite angular cell bound."); + // Both signs share the derivative bounds; only the center value changes. + upper[c] = std::max(0.0, center + excursion); + negative_upper[c] = std::max(0.0, -center + excursion); + } + } + +private: + int nt_, nz_, m_, n_; + vector theta_, zeta_, inverse_; +}; + // Bernstein basis functions are nonnegative and sum to one on [0,1], so the // largest Bernstein coefficient bounds the cubic for rejection sampling. double cubic_rejection_bound(const array& c) @@ -1330,8 +1413,8 @@ void StellaratorSource::precompute_sampling_distributions() // of poloidal degree <= 3*m_max and toroidal degree <= 3*n_max per field // period, so a uniform trapezoidal product rule with more points than the // bandwidth integrates it exactly (discrete Fourier orthogonality). The - // denser grid below also estimates the angular maximum for rejection - // sampling; a heuristic safety factor allows for peaks between grid points. + // denser grid below also resolves every density Fourier coefficient used + // to construct the angular cell bounds (more than twice the bandwidth). int m_max = 0; int n_max = 0; for (size_t k = 0; k < n_modes; ++k) { @@ -1414,8 +1497,15 @@ void StellaratorSource::precompute_sampling_distributions() double quad_wgt = 4.0 * PI * PI / n_grid; // full-torus angular weight double v_tot = 0.0; vector vp_lo(n_bins), vp_mid(n_bins), vp_hi(n_bins); - vector negative_envelope(n_bins); - envelope_.assign(n_bins, 0.0); + constexpr size_t n_cells = ANGULAR_BINS * ANGULAR_BINS; + angular_bounds_.resize({n_bins, size_t {3}, n_cells}); + tensor::Tensor negative_bounds({n_bins, size_t {3}, n_cells}); + angular_cells_.resize(n_bins); + angular_totals_.resize(n_bins); + AngularBounds bounds(n_theta, n_zeta, 3 * m_max, 3 * n_max); + array, 3> controls; + for (auto& v : controls) + v.resize(n_grid); double f_peak = 0.0; double f_min = 0.0; evaluate_surface(0); @@ -1433,8 +1523,6 @@ void StellaratorSource::precompute_sampling_distributions() const double* Zt0 = &Zt_s[lo]; const double* Zt1 = &Zt_s[hi]; double s0 = 0.0, sh = 0.0, s1 = 0.0; - double env = 0.0; - double negative_env = 0.0; for (size_t p = 0; p < n_grid; ++p) { // Linear interpolation makes radial derivatives constant within a bin. // R and the poloidal derivatives are linear, so R*tau is quadratic. @@ -1449,11 +1537,11 @@ void StellaratorSource::precompute_sampling_distributions() s0 += f0; sh += fh; s1 += f1; - - // Endpoints and midpoint determine the radial quadratic at this angle. - // Bound both signs until the global Jacobian orientation is known. - env = std::max(env, quad_max01(f0, fh, f1)); - negative_env = std::max(negative_env, quad_max01(-f0, -fh, -f1)); + // Express the radial quadratic in the Bernstein basis so its cell + // bounds can be mixed with nonnegative weights at the sampled radius. + controls[0][p] = f0; + controls[1][p] = 2 * fh - (f0 + f1) / 2; + controls[2][p] = f1; f_peak = std::max({f_peak, f0, fh, f1}); f_min = std::min({f_min, f0, fh, f1}); } @@ -1461,14 +1549,15 @@ void StellaratorSource::precompute_sampling_distributions() vp_mid[b] = sh * quad_wgt; vp_hi[b] = s1 * quad_wgt; v_tot += sh * quad_wgt * (rho_[b + 1] - rho_[b]); - // Heuristic allowance for angular peaks between quadrature points - envelope_[b] = 1.05 * env; - negative_envelope[b] = 1.05 * negative_env; + for (int j = 0; j < 3; ++j) { + bounds.compute( + controls[j], &angular_bounds_(b, j, 0), &negative_bounds(b, j, 0)); + } } jacobian_sign_ = (v_tot >= 0.0) ? 1.0 : -1.0; if (jacobian_sign_ < 0.0) { // Reverse the integrals and extrema to use a positive volume element. - envelope_.swap(negative_envelope); + angular_bounds_ = std::move(negative_bounds); double old_peak = f_peak; f_peak = -f_min; f_min = -old_peak; @@ -1485,6 +1574,19 @@ void StellaratorSource::precompute_sampling_distributions() "coefficient tables are inconsistent."); } + // Build cell-selection tables after choosing the global orientation. + for (size_t b = 0; b < n_bins; ++b) { + for (int j = 0; j < 3; ++j) { + double* upper = &angular_bounds_(b, j, 0); + double total = std::accumulate(upper, upper + n_cells, 0.0); + if (!std::isfinite(total)) + fatal_error("StellaratorSource: nonfinite angular bound integral."); + angular_totals_[b][j] = total; + // An identically zero component (e.g., on the axis) is never selected. + if (total > 0) + angular_cells_[b][j].assign({upper, n_cells}); + } + } // Integrate S(rho)*V'(rho) independently on each radial interval. Scaling // the emission profile by its maximum preserves the normalized distribution. double s_max = @@ -1639,13 +1741,35 @@ SourceSite StellaratorSource::sample(uint64_t* seed) const double rho = rho_[bin] + t * (rho_[bin + 1] - rho_[bin]); // 3. Rejection-sample (theta, zeta) from p(theta, zeta | rho) ~ R*tau - double env = envelope_[bin]; + array beta {(1 - t) * (1 - t), 2 * t * (1 - t), t * t}; + array weights; + for (int j = 0; j < 3; ++j) + weights[j] = beta[j] * angular_totals_[bin][j]; + double total = weights[0] + weights[1] + weights[2]; + if (!(total > 0) || !std::isfinite(total)) + fatal_error("StellaratorSource: invalid angular proposal normalization."); double R, Z, zeta; n_reject = 0; while (true) { - double theta = uniform_distribution(0.0, 2.0 * PI, seed); - zeta = uniform_distribution(0.0, 2.0 * PI, seed); + // Mix the three alias tables to select cells proportional to their + // quadratic bound at this radius. Resample the cell on every rejection. + double u = prn(seed) * total; + int j = u < weights[0] ? 0 : (u < weights[0] + weights[1] ? 1 : 2); + size_t cell = angular_cells_[bin][j].sample(seed); + double theta = (cell / ANGULAR_BINS + prn(seed)) * (2 * PI / ANGULAR_BINS); + zeta = (cell % ANGULAR_BINS + prn(seed)) * + (2 * PI / (ANGULAR_BINS * num_field_periods_)); + double env = beta[0] * angular_bounds_(bin, 0, cell) + + beta[1] * angular_bounds_(bin, 1, cell) + + beta[2] * angular_bounds_(bin, 2, cell); double f = eval_density(bin, t, theta, zeta, &R, &Z); + if (!std::isfinite(f) || f < 0) + fatal_error("StellaratorSource: invalid angular density; check the " + "Fourier geometry."); + if (f > env) + fatal_error(fmt::format("StellaratorSource: angular density {} exceeds " + "cell bound {} in radial interval {}.", + f, env, bin)); if (prn(seed) * env < f) break; if (++n_reject > MAX_SOURCE_REJECTIONS_PER_SAMPLE) { @@ -1655,6 +1779,8 @@ SourceSite StellaratorSource::sample(uint64_t* seed) const } // 4. Convert to Cartesian coordinates (zeta is the cylindrical angle) + zeta += 2 * PI * static_cast(num_field_periods_ * prn(seed)) / + num_field_periods_; site.r = {R * std::cos(zeta), R * std::sin(zeta), Z}; // 5. Sample isotropic direction diff --git a/tests/unit_tests/test_source_stellarator.py b/tests/unit_tests/test_source_stellarator.py index 7e8e28b0859..f18fe4197f9 100644 --- a/tests/unit_tests/test_source_stellarator.py +++ b/tests/unit_tests/test_source_stellarator.py @@ -194,6 +194,48 @@ def test_stellarator_source_axisymmetric_sampling(run_in_tmpdir): assert_sample_mean(major_r, r0 + a**2 / (4.0 * r0)) +@pytest.mark.parametrize('nfp, orientation', [(1, 1), (3, -1)]) +def test_stellarator_source_angular_cells(run_in_tmpdir, nfp, orientation): + """A phase-shifted tube has analytic angular and radial marginals. + + Its narrow peaks exceed the old grid maximum plus 5% envelope. Two radial + surfaces exercise the zero axis component and the full quadratic mixture; + toroidal density harmonics exceed the cell grid's Nyquist frequency. + """ + r0, a, epsilon, k, phase = 600., 40., .9, 8, np.pi / 8 + rho = np.array([0., 1.]) + rc = np.array([[r0, 0., 0., 0.], + [r0, a, a * epsilon * np.cos(phase) / 2, + a * epsilon * np.cos(phase) / 2]]) + zs = rc.copy() + zs[:, 0] = 0. + rs = np.zeros_like(rc) + rs[1, 2:] = a * epsilon * np.sin(phase) / 2 * np.array([-1, 1]) + zc = -rs + src = make_source( + rho=rho, emission_density=np.ones(2), mode_m=[0, 1, 1, 1], + mode_n=[0, 0, k, -k], rmnc=rc, zmns=orientation * zs, + rmns=rs, zmnc=orientation * zc, num_field_periods=nfp, + energy=openmc.stats.delta_function(14.07e6), + ) + sites = make_model(src).sample_external_source(100_000) + xyz = np.array([s.r for s in sites]) + r = np.hypot(xyz[:, 0], xyz[:, 1]) + phi = np.arctan2(xyz[:, 1], xyz[:, 0]) + harmonic = np.cos(k * nfp * phi - phase) + minor_rho = np.hypot(r - r0, xyz[:, 2]) / (a * (1 + epsilon * harmonic)) + assert minor_rho.max() <= 1 + 1e-10 + assert_sample_mean(minor_rho, 2 / 3) + assert_sample_mean(harmonic, epsilon / (1 + epsilon**2 / 2)) + assert_sample_mean(np.sin(k * nfp * phi - phase), 0.) + assert_sample_mean(r, r0 + a**2 * (1 + 3 * epsilon**2 + + 3 * epsilon**4 / 8) / (4 * r0 * (1 + epsilon**2 / 2))) + if nfp > 1: + # Sampling only one field period would fail these full-torus moments. + assert_sample_mean(np.cos(phi), 0.) + assert_sample_mean(np.sin(phi), 0.) + + @pytest.mark.parametrize('orientation', [1, -1]) def test_stellarator_source_rotating_ellipse_sampling(run_in_tmpdir, orientation): """Compare sampled moments of a 3-D configuration against quadrature of From 6069bee79146457666d745a02540cfcdc094e47b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 10 Oct 2026 14:11:06 -0500 Subject: [PATCH 02/15] Use a reproducible seed in angular cell sampling tests --- tests/unit_tests/test_source_stellarator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_source_stellarator.py b/tests/unit_tests/test_source_stellarator.py index f18fe4197f9..9705436add5 100644 --- a/tests/unit_tests/test_source_stellarator.py +++ b/tests/unit_tests/test_source_stellarator.py @@ -218,7 +218,7 @@ def test_stellarator_source_angular_cells(run_in_tmpdir, nfp, orientation): rmns=rs, zmnc=orientation * zc, num_field_periods=nfp, energy=openmc.stats.delta_function(14.07e6), ) - sites = make_model(src).sample_external_source(100_000) + sites = make_model(src).sample_external_source(100_000, prn_seed=20261010) xyz = np.array([s.r for s in sites]) r = np.hypot(xyz[:, 0], xyz[:, 1]) phi = np.arctan2(xyz[:, 1], xyz[:, 0]) From 37a0b2cd72fd79660cd67b0e3caab0185d94d6c4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Oct 2026 11:50:34 -0500 Subject: [PATCH 03/15] Add source sampling docs --- docs/source/methods/index.rst | 1 + docs/source/methods/source_sampling.rst | 720 ++++++++++++++++++++++++ docs/source/usersguide/settings.rst | 9 + 3 files changed, 730 insertions(+) create mode 100644 docs/source/methods/source_sampling.rst diff --git a/docs/source/methods/index.rst b/docs/source/methods/index.rst index 121d04b1ded..071ee2dc7ed 100644 --- a/docs/source/methods/index.rst +++ b/docs/source/methods/index.rst @@ -12,6 +12,7 @@ Theory and Methodology geometry cross_sections random_numbers + source_sampling neutron_physics photon_physics charged_particles_physics diff --git a/docs/source/methods/source_sampling.rst b/docs/source/methods/source_sampling.rst new file mode 100644 index 00000000000..fa92a4f13a8 --- /dev/null +++ b/docs/source/methods/source_sampling.rst @@ -0,0 +1,720 @@ +.. _methods_source_sampling: + +=============== +Source Sampling +=============== + +A source site specifies a particle's position :math:`\mathbf{r}`, direction +:math:`\mathbf{\Omega}`, energy :math:`E`, emission time :math:`t`, particle +type, and statistical weight :math:`w`. In a fixed source calculation, OpenMC +draws sites from the external source throughout the calculation. In an +eigenvalue calculation, the external source supplies the initial population; +subsequent generations use the fission bank as described in +:ref:`method-successive-generations`. + +This chapter describes the algorithms represented by the subclasses of +:class:`openmc.SourceBase`. Input examples are given in +:ref:`usersguide_source`. Unless otherwise stated, :math:`\xi` denotes an +independent uniform random number on :math:`[0,1)`, generated as described in +:ref:`methods_random_numbers`. Probability densities below are normalized with +respect to the coordinates explicitly indicated. + +.. _methods_source_selection: + +---------------------------- +Selecting Among Source Types +---------------------------- + +Suppose there are :math:`N` external sources with nonnegative strengths +:math:`s_i`, with :math:`\sum_i s_i>0`, and normalized distributions +:math:`p_i`. Their mixture is + +.. math:: + :label: source-mixture + + p(\mathbf{r},\mathbf{\Omega},E,t) + = \sum_{i=1}^{N} P_i p_i(\mathbf{r},\mathbf{\Omega},E,t), + \qquad P_i = \frac{s_i}{\sum_j s_j}. + +OpenMC first selects a source index with probability :math:`P_i` and then calls +that source's sampling method. Discrete selections of this kind use alias +tables: after initialization, each selection requires a fixed number of random +draws and table lookups, independent of the number of entries. + +With :attr:`~openmc.Settings.uniform_source_sampling` enabled, the source index +is instead selected with probability :math:`1/N`. The sampled site's weight is +multiplied by + +.. math:: + + \frac{P_i}{1/N} = \frac{N s_i}{\sum_j s_j}. + +This preserves the expected contribution of each source while allocating the +same expected number of histories to every source. These factors multiply any +weights supplied by the selected source. Source strengths determine the +relative probabilities; converting results per source particle to an absolute +physical rate is a separate normalization described in +:ref:`usersguide_tally_normalization`. + +------------------------- +Independent Distributions +------------------------- + +An :class:`openmc.IndependentSource` defines separate distributions for +position, direction, energy, and time. Before applying constraints, its joint +distribution factors as + +.. math:: + + p(\mathbf{r},\mathbf{\Omega},E,t) + = p_r(\mathbf{r}) p_\Omega(\mathbf{\Omega}) p_E(E) p_t(t). + +The particle type is fixed for each independent source. Each factor can itself +be continuous, discrete, or a mixture of distributions. + +Sampling One-Dimensional Distributions +-------------------------------------- + +For a normalized density :math:`f(x)`, inverse-transform sampling draws +:math:`x=F^{-1}(\xi)`, where + +.. math:: + + F(x) = \int_{-\infty}^{x} f(u)\,du. + +OpenMC provides analytic distributions as well as tabulated distributions +through :mod:`openmc.stats`. A tabulated distribution first selects an interval +using its cumulative probability, then inverts the CDF within that interval +according to the specified interpolation. For example, a linear-linear PDF on +:math:`[x_i,x_{i+1}]` has + +.. math:: + + f(x_i+y) = f_i + m_i y, \qquad + F(x_i+y)-F(x_i) = f_i y + \frac{m_i}{2}y^2, + +where :math:`m_i=(f_{i+1}-f_i)/(x_{i+1}-x_i)`. Inverting this quadratic gives +the sampled value. A histogram PDF has a linear CDF within each interval. +Discrete distributions select among their supplied values using their +probabilities; mixture distributions select a component and sample it. + +Energy and time distributions use the same one-dimensional sampling +infrastructure. Analytic spectra, such as Maxwell and Watt spectra, use their +specialized sampling algorithms. A :class:`~openmc.stats.DecaySpectrum` selects +from the nuclide decay spectra supplied by the depletion chain, sets the source +strength from the total photon emission rate, and records the parent nuclide +with each sampled photon. + +Spatial Coordinates and Volume Factors +-------------------------------------- + +A :class:`~openmc.stats.CartesianIndependent` distribution samples :math:`x`, +:math:`y`, and :math:`z` independently. A +:class:`~openmc.stats.CylindricalIndependent` distribution samples +:math:`(R,\phi,z)` and transforms to Cartesian coordinates. Its supplied radial +PDF is a density with respect to :math:`dR`. For a uniform density per unit +volume in a cylindrical shell, the volume element +:math:`dV=R\,dR\,d\phi\,dz` requires + +.. math:: + + p_R(R) = \frac{2R}{R_{\max}^2-R_{\min}^2}, \qquad + R = \sqrt{R_{\min}^2 + \xi(R_{\max}^2-R_{\min}^2)}. + +Similarly, :class:`~openmc.stats.SphericalIndependent` samples radius, +:math:`\mu=\cos\theta`, and azimuth. Uniform volume sampling in a spherical +shell requires :math:`p_r(r)\propto r^2`, giving + +.. math:: + + r = \left[r_{\min}^3 + \xi(r_{\max}^3-r_{\min}^3)\right]^{1/3}. + +The coordinate distributions specify these radial factors explicitly. The +coordinate transformation itself introduces no additional sampling weight. +Points and boxes are special cases with fixed coordinates or uniform Cartesian +coordinates, respectively. + +A :class:`~openmc.stats.PointCloud` is a discrete spatial distribution. Given +coordinates :math:`\mathbf{r}_j` with strengths :math:`a_j`, it selects +:math:`\mathbf{r}_j` with probability :math:`a_j/\sum_k a_k`. It places the +particle exactly at the selected point. If points approximate a continuous +emission density, their strengths must incorporate the volumes represented by +the points. + +Angular Distributions +--------------------- + +An isotropic direction has density :math:`1/(4\pi)` per unit solid angle. +Because :math:`d\Omega=d\mu\,d\varphi`, it can be sampled using + +.. math:: + + \mu = 2\xi_1-1, \qquad \varphi=2\pi\xi_2, + +and constructing the direction from the polar cosine and azimuth. A +:class:`~openmc.stats.PolarAzimuthal` distribution supplies separate PDFs for +:math:`\mu` and :math:`\varphi`, referred to a specified orientation, while +:class:`~openmc.stats.Monodirectional` always returns one direction. These +angles describe particle flight directions and are distinct from the spatial +angles used to locate a source site. + +------------------------- +Sources Defined on Meshes +------------------------- + +A mesh distribution first selects an element, then draws a position uniformly +within that element. If :math:`a_j` is the integrated strength assigned to +element :math:`j`, with volume :math:`V_j`, the spatial PDF is + +.. math:: + + p_r(\mathbf{r}) = \frac{a_j}{\sum_k a_k}\frac{1}{V_j}, + \qquad \mathbf{r}\in V_j. + +An :class:`openmc.IndependentSource` with a +:class:`~openmc.stats.MeshSpatial` distribution uses the same direction, energy, +and time distributions for all elements. With ``volume_normalized=True``, the +supplied strengths represent densities per unit volume: OpenMC multiplies them +by element volumes before selecting an element. Otherwise, the supplied values +are integrated element strengths. + +An :class:`openmc.MeshSource` assigns an independent source to each element. +It selects an element using that source's strength, replaces the source's +spatial distribution with uniform sampling inside the selected element, and +then samples the element's direction, energy, time, and particle type. This +allows those distributions to vary with position through the element index. +Its total strength is the sum of the element source strengths. Element +strengths are integrated strengths; a prescribed volumetric emission density +must be multiplied by element volume when constructing them. + +--------------------------- +Stored and Compiled Sources +--------------------------- + +An :class:`openmc.FileSource` loads source sites from a supported particle file. +If the file contains :math:`M` sites, each draw selects an index uniformly from +:math:`0,\ldots,M-1` and copies that site's position, direction, energy, time, +weight, and particle type. Sampling is with replacement, so repeated sites are +possible. Correlations stored in a site, such as those between position and +energy, are preserved. The stored weight is retained rather than used as the +probability of choosing that site. + +An :class:`openmc.CompiledSource` calls a sampling method in a user-supplied +shared library. That method returns a complete source site and can implement +correlations or distributions beyond those exposed by the built-in classes. +Its sampling and statistical-weight corrections are the responsibility of the +custom implementation. Source selection and applicable source constraints +surround the custom sampling call as they do for other source types. + +.. _methods_plasma_sources: + +--------------------- +Plasma Volume Sources +--------------------- + +The :class:`openmc.TokamakSource` and :class:`openmc.StellaratorSource` classes +sample neutron emission from flux surfaces. Both assume a nonnegative emission +density :math:`S` per unit physical volume that depends only on a radial flux +coordinate. This is an assumption of these source models; an equilibrium's +geometry alone does not determine its neutron emission profile. + +For coordinates :math:`(u,\alpha,\phi)` with volume Jacobian +:math:`J=|\partial(x,y,z)/\partial(u,\alpha,\phi)|`, the spatial density in +those coordinates is + +.. math:: + :label: plasma-source-joint + + p(u,\alpha,\phi) \propto S(u)J(u,\alpha,\phi). + +Thus uniform angular coordinates generally do not give uniform emission per +unit volume. Define the differential volume + +.. math:: + + V'(u) = \int\!\int J(u,\alpha,\phi)\,d\alpha\,d\phi. + +The joint PDF factors into a radial marginal and an angular conditional: + +.. math:: + :label: plasma-source-factorization + + p_u(u) = \frac{S(u)V'(u)}{\int S(v)V'(v)\,dv}, \qquad + p(\alpha,\phi\mid u) = \frac{J(u,\alpha,\phi)}{V'(u)}. + +Sampling in this order accounts for the volume occupied by each flux interval. +The emission profile cancels from the conditional angular PDF. + +Radially Varying Spectra +------------------------ + +Both plasma sources accept either one energy distribution or an energy +distribution at every radial grid point. In an interval with local coordinate +:math:`\eta\in[0,1]`, the conditional energy distribution is the mixture + +.. math:: + + p_E(E\mid u) = (1-\eta)p_{E,i}(E) + \eta p_{E,i+1}(E). + +OpenMC selects distribution :math:`i+1` with probability :math:`\eta` and +distribution :math:`i` otherwise, then samples energy from the selected +distribution. This interpolates the distributions, rather than interpolating +two sampled energies. Both classes sample an isotropic flight direction and a +separate time distribution, which defaults to a point mass at zero. Importance +weights returned by the energy and time distributions multiply the site's +weight. + +.. _methods_tokamak_source: + +Tokamak Sources +--------------- + +The tokamak geometry uses the parameterization + +.. math:: + + \begin{aligned} + R(r,\alpha) &= R_0+r\cos\psi+\Delta[1-(r/a)^2], \\ + Z(r,\alpha) &= Z_{\mathrm{shift}}+\kappa r\sin\alpha, \\ + \psi &= \alpha+\delta\sin\alpha, + \end{aligned} + +where :math:`a` is the minor radius, :math:`\kappa` is elongation, +:math:`\delta` controls triangularity, and :math:`\Delta` is the Shafranov +shift. The cylindrical azimuth :math:`\phi` is independent of these functions. +Let :math:`q=r/a`, :math:`\epsilon=a/R_0`, and :math:`D=\Delta/a`. +Differentiating the geometry gives + +.. math:: + + \begin{aligned} + R_r Z_\alpha-R_\alpha Z_r + &= \kappa r\,[G(\alpha)-2Dq\cos\alpha], \\ + G(\alpha) + &= \cos(\delta\sin\alpha) + +\delta\sin\alpha\cos\alpha\sin\psi. + \end{aligned} + +The coordinate volume element is :math:`R` times this poloidal Jacobian. +Up to constants independent of :math:`q` and :math:`\alpha`, the joint density +is therefore + +.. math:: + :label: tokamak-source-joint + + p(q,\alpha) \propto S(q)q + [1+\epsilon q\cos\psi+\epsilon D(1-q^2)] + [G(\alpha)-2Dq\cos\alpha]. + +Sampling assumes a usable, consistently oriented map. The elongation +contributes a constant to the volume Jacobian and cancels from the normalized +PDF. The toroidal angle is uniform over the requested interval. + +Radial Marginal +~~~~~~~~~~~~~~~ + +Integrating :eq:`tokamak-source-joint` over the poloidal angle yields + +.. math:: + :label: tokamak-source-radial + + p_q(q) \propto S(q)q H(q), \qquad + H(q) = 1+\epsilon D-\frac{3}{8}c_1\epsilon q-2\epsilon Dq^2, + +where, in terms of Bessel functions of the first kind, + +.. math:: + + c_0 = J_0(\delta)+J_2(\delta), \qquad + c_1 = \frac{J_1(2\delta)+J_3(2\delta)}{c_0}. + +The emission profile is linearly interpolated on the user grid. OpenMC refines +each interval into at least eight subintervals with spacing at most +:math:`10^{-3}`, evaluates :eq:`tokamak-source-radial` there, and constructs a +linear-linear tabulated distribution. Radial sampling then uses CDF inversion +of that tabulation. For circular, concentric surfaces, +:math:`\delta=\Delta=0`, the radial marginal reduces to +:math:`p_q(q)\propto S(q)q`. + +Conditional Poloidal Sampling +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The two bracketed factors in :eq:`tokamak-source-joint` are quadratic and +linear in :math:`q`. Write their Bernstein forms as + +.. math:: + + \begin{aligned} + \frac{R}{R_0} &= b_0(1-q)^2+2b_1q(1-q)+b_2q^2, \\ + G-2Dq\cos\alpha &= b_3(1-q)+b_4q, + \end{aligned} + +with + +.. math:: + + \begin{aligned} + b_0 &= 1+\epsilon D, & b_1 &= b_0+\tfrac12\epsilon\cos\psi, \\ + b_2 &= 1+\epsilon\cos\psi, & b_3 &= G(\alpha), \\ + b_4 &= G(\alpha)-2D\cos\alpha. + \end{aligned} + +Their product is a sum of six separable terms +:math:`\sum_{k=0}^5 w_k(q)g_k(\alpha)`, where + +.. math:: + + \begin{aligned} + (g_0,g_1,g_2,g_3,g_4,g_5) + &= (b_0b_3,b_1b_3,b_2b_3,b_0b_4,b_1b_4,b_2b_4), \\ + (w_0,w_1,w_2,w_3,w_4,w_5) + &= ((1-q)^3,\,2q(1-q)^2,\,q^2(1-q), \\ + &\qquad q(1-q)^2,\,2q^2(1-q),\,q^3). + \end{aligned} + +Define :math:`I_k=\int_0^\pi g_k(\alpha)\,d\alpha` and scaled integrals +:math:`\widehat I_k=I_k/(\pi c_0)`. Analytic integration gives + +.. math:: + + \begin{aligned} + \widehat I_0 &= 1+\epsilon D, \\ + \widehat I_1 &= 1+\epsilon D-\tfrac{3}{16}c_1\epsilon, \\ + \widehat I_2 &= 1-\tfrac38 c_1\epsilon, \\ + \widehat I_3 &= 1+\epsilon D, \\ + \widehat I_4 &= 1+\tfrac12\epsilon D-\tfrac{3}{16}c_1\epsilon, \\ + \widehat I_5 &= 1-\epsilon D-\tfrac38 c_1\epsilon. + \end{aligned} + +Consequently, the conditional PDF on :math:`[0,\pi]` is a mixture of the six +normalized densities :math:`g_k/I_k`, with component probabilities +:math:`w_k(q)I_k/\sum_j w_j(q)I_j`. The implementation tabulates each component +using ``n_alpha`` points and samples its tabulated CDF after component +selection. Up-down symmetry supplies the other half-plane by +replacing :math:`\alpha` with :math:`2\pi-\alpha` with probability one half. +The position is then mapped to :math:`(R\cos\phi,R\sin\phi,Z)`. This +construction samples the spatial coordinates without a rejection loop, with +discretization error from the radial and poloidal tabulations. + +.. note:: + + The continuous integrals above give the required mixture weights. The + current implementation stores different values for :math:`I_2`, + :math:`I_4`, and :math:`I_5`; that mismatch introduces an additional error + in poloidal sampling beyond the tabulation error. + +.. _methods_stellarator_source: + +Stellarator Sources +------------------- + +Fourier Representation and Volume Jacobian +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The stellarator source uses flux coordinates :math:`(\rho,\theta,\zeta)`, +where :math:`\rho=\sqrt{s}` for normalized toroidal flux :math:`s`, +:math:`\theta` is a poloidal angle, and :math:`\zeta` equals the cylindrical +azimuth. With :math:`N_{\mathrm{fp}}` field periods, define +:math:`\chi_k=m_k\theta-n_kN_{\mathrm{fp}}\zeta`. The surfaces are + +.. math:: + :label: stellarator-source-fourier + + \begin{aligned} + R(\rho,\theta,\zeta) + &= \sum_k [R_k^c(\rho)\cos\chi_k+R_k^s(\rho)\sin\chi_k], \\ + Z(\rho,\theta,\zeta) + &= \sum_k [Z_k^s(\rho)\sin\chi_k+Z_k^c(\rho)\cos\chi_k]. + \end{aligned} + +The sine coefficients of :math:`R` and cosine coefficients of :math:`Z` vanish +for stellarator-symmetric geometry. The same sampling procedure handles both +symmetric and asymmetric tables. + +The Cartesian position is +:math:`\mathbf{r}=(R\cos\zeta,R\sin\zeta,Z)`. Applying the chain rule through +cylindrical coordinates gives + +.. math:: + :label: stellarator-source-jacobian + + \left|\frac{\partial(x,y,z)}{\partial(\rho,\theta,\zeta)}\right| + = R|\tau|, \qquad + \tau=R_\rho Z_\theta-R_\theta Z_\rho. + +Here subscripts denote partial derivatives. The derivatives :math:`R_\zeta` +and :math:`Z_\zeta` drop out of the determinant because the cylindrical azimuth +is exactly :math:`\zeta`. Thus the joint coordinate PDF is +:math:`p(\rho,\theta,\zeta)\propto S(\rho)R|\tau|`. + +The map is assumed to have positive :math:`R`, nested surfaces, and a +consistently oriented Jacobian away from the magnetic axis. Either orientation +is allowed. Initialization determines a global sign :math:`\sigma` from a +signed volume estimate and evaluates :math:`f=\sigma R\tau`, so that +:math:`f=R|\tau|` for a valid map. Grid checks reject nonpositive radii or +significant changes in Jacobian sign; these checks are not a global geometric +certification. + +Interpolation Within a Radial Interval +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +On interval :math:`i`, define its width :math:`h_i=\rho_{i+1}-\rho_i` and local +coordinate :math:`\eta=(\rho-\rho_i)/h_i`. Each Fourier coefficient :math:`C` +and the emission profile are interpolated linearly: + +.. math:: + + C(\eta)=(1-\eta)C_i+\eta C_{i+1}, \qquad + S(\eta)=S_i+\eta(S_{i+1}-S_i). + +Because Fourier evaluation is linear in the coefficients, at fixed angles +:math:`R`, :math:`Z`, :math:`R_\theta`, and :math:`Z_\theta` are linear in +:math:`\eta`, while + +.. math:: + + R_\rho=\frac{R_{i+1}-R_i}{h_i}, \qquad + Z_\rho=\frac{Z_{i+1}-Z_i}{h_i} + +are constant within the interval. It follows that :math:`\tau` is linear in +:math:`\eta` and :math:`f=\sigma R\tau` is quadratic. Angular integration +preserves this polynomial degree, so the differential volume within interval +:math:`i` is also quadratic: + +.. math:: + :label: stellarator-source-volume-polynomial + + V'_i(\eta)=A_0+A_1\eta+A_2\eta^2. + +The derivative generally has different one-sided values at an interior radial +surface because the coefficient slopes change there. Each interval uses its +own values of :math:`R_\rho` and :math:`Z_\rho`; the endpoint volumes from +adjacent intervals are kept separate. + +Computing Differential Volumes +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Initialization evaluates :math:`f` at :math:`\eta=0,1/2,1` on a uniform +angular grid. Its sizes are + +.. math:: + + N_\theta=\max(64,8m_{\max}), \qquad + N_\zeta=\max(64,8|n|_{\max}). + +The poloidal grid covers :math:`[0,2\pi)` and the toroidal grid covers one field +period, :math:`[0,2\pi/N_{\mathrm{fp}})`. The full-torus integral is + +.. math:: + + V'_i(\eta) + = \frac{4\pi^2}{N_\theta N_\zeta} + \sum_{j=0}^{N_\theta-1}\sum_{k=0}^{N_\zeta-1} + f_i(\eta,\theta_j,\zeta_k). + +The weight includes all field periods. Products in :math:`R\tau` have poloidal +and toroidal degrees at most three times the respective geometry mode limits. +The trapezoidal grid therefore integrates this finite Fourier series exactly +in exact arithmetic, provided the Jacobian has a consistent sign. + +Let :math:`v_0`, :math:`v_m`, and :math:`v_1` be these integrals at +:math:`\eta=0,1/2,1`. Interpolation of the quadratic gives + +.. math:: + + A_0=v_0, \qquad + A_1=4v_m-3v_0-v_1, \qquad + A_2=2(v_0+v_1-2v_m). + +Only two adjacent surfaces' angular grids are held in memory at once. Their +positions and poloidal derivatives give the interval's radial derivatives by +finite differences. Initialization retains the per-interval integrals and +angular bounds needed for sampling. + +Selecting and Sampling a Radial Interval +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Multiplying the linear emission profile by +:eq:`stellarator-source-volume-polynomial` yields an unnormalized cubic radial +PDF in the local coordinate: + +.. math:: + :label: stellarator-source-radial-polynomial + + \begin{aligned} + P_i(\eta) &= S(\eta)V'_i(\eta) + = c_0+c_1\eta+c_2\eta^2+c_3\eta^3, \\ + c_0 &= S_i A_0, \\ + c_1 &= S_i A_1+\Delta S_i A_0, \\ + c_2 &= S_i A_2+\Delta S_i A_1, \\ + c_3 &= \Delta S_i A_2, + \end{aligned} + +where :math:`\Delta S_i=S_{i+1}-S_i`. Its integrated mass is + +.. math:: + :label: stellarator-source-interval-mass + + \begin{aligned} + M_i &= h_i\int_0^1 P_i(\eta)\,d\eta \\ + &= h_i\left(\frac{S_i v_0}{6} + +\frac{(S_i+S_{i+1})v_m}{3} + +\frac{S_{i+1}v_1}{6}\right). + \end{aligned} + +The last expression is Simpson's rule, which is exact for a cubic. OpenMC +selects interval :math:`i` with probability :math:`M_i/\sum_j M_j` using an +alias table. A common scaling of all emission values cancels from these +probabilities; the implementation divides them by their maximum before forming +the masses and PDFs. + +Within the selected interval, OpenMC uses uniform proposals +:math:`\eta=\xi_1` with rejection sampling. An upper bound follows by writing +the cubic in the Bernstein basis: + +.. math:: + + P_i(\eta) = \sum_{j=0}^{3}\beta_j + {3\choose j}\eta^j(1-\eta)^{3-j}, + +where + +.. math:: + + \begin{aligned} + \beta_0 &= c_0, \\ + \beta_1 &= c_0+c_1/3, \\ + \beta_2 &= c_0+2c_1/3+c_2/3, \\ + \beta_3 &= c_0+c_1+c_2+c_3. + \end{aligned} + +The Bernstein basis functions are nonnegative and sum to one, so +:math:`B_i^{r}=\max_j\beta_j` bounds :math:`P_i` throughout the interval. +The implementation adds a small relative margin for floating-point roundoff +and stores the coefficients divided by that bound. A proposal is accepted when + +.. math:: + + \xi_2 < \frac{P_i(\eta)}{B_i^{r}}. + +On rejection, only :math:`\eta` is redrawn; the selected interval is retained. +Selecting a new interval after a rejection would change its probability in +proportion to its acceptance rate. After acceptance, +:math:`\rho=\rho_i+h_i\eta` has the required marginal +:math:`p_\rho\propto S V'`. + +Conditional Angular Rejection +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +With :math:`\rho` fixed, the target PDF is + +.. math:: + + p(\theta,\zeta\mid\rho) + = \frac{f_i(\eta,\theta,\zeta)}{V'_i(\eta)}. + +At each angular grid point, the three radial evaluations determine the +quadratic :math:`f_i(\eta)` throughout the interval. Its maximum is obtained +from the endpoints and, for a concave quadratic, an interior vertex. The current +angular envelope is + +.. math:: + + B_i^{\mathrm{ang}} = 1.05\max_{j,k}\max_{0\le\eta\le1} + f_i(\eta,\theta_j,\zeta_k). + +The radial maximization here is exact in exact arithmetic. The angular maximum +is estimated on a grid, with a 5% heuristic allowance for peaks between nodes. +Exact Fourier quadrature of an integral does not establish an upper bound on +the function. This envelope therefore does not provide a rigorous bound for +every supported coefficient table; correctness of angular rejection requires +it to bound the continuous angular density. + +The sampler proposes :math:`\theta=2\pi\xi_1` and +:math:`\zeta=2\pi\xi_2` over the full torus, and accepts when + +.. math:: + + \xi_3 B_i^{\mathrm{ang}} < f_i(\eta,\theta,\zeta). + +Rejected angles are redrawn while retaining both the radial interval and +:math:`\rho`. If the envelope is valid, this produces the angular conditional +PDF above. An underestimated envelope clips the peaks because acceptance +cannot exceed one; the current sampling loop has no explicit check for that +condition. The angular rejection limit detects excessive rejection, which is a +different failure condition. + +Sampling Sequence and Approximation +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The steps in the stellarator source's sampling method are: + +1. Select a radial interval using its integrated mass :math:`M_i`. +2. Sample the local radial coordinate by cubic-PDF rejection within that + interval, and calculate :math:`\rho`. +3. Sample the poloidal and toroidal angles by conditional rejection at that + fixed radius. +4. Map the accepted position to :math:`(R\cos\zeta,R\sin\zeta,Z)`. +5. Sample an isotropic flight direction. +6. Sample energy from the common spectrum or the mixture of spectra bracketing + :math:`\rho`, then sample emission time and apply their weight factors. + +Radial integration and rejection sampling are consistent with the piecewise +linear coefficient and emission tables, up to floating-point effects. Those +tables approximate the underlying equilibrium and emission profile; increasing +the number of radial surfaces can reduce interpolation error. The finite +Fourier representation also limits the geometric detail retained. These input +approximations are distinct from the angular-envelope condition above. + +The :meth:`~openmc.StellaratorSource.from_vmec` and +:meth:`~openmc.StellaratorSource.from_desc` constructors prepare the coefficient +tables before C++ initialization. VMEC's surfaces, tabulated in :math:`s`, are +relabeled by :math:`\rho=\sqrt{s}`. For DESC, the Fourier-Zernike radial basis +is evaluated at the requested radial surfaces and the angular product basis is +converted to :eq:`stellarator-source-fourier`. Both constructors convert +geometry coefficients from meters to centimeters. The user supplies the +emission profile and energy distributions separately. + +.. _methods_source_constraints: + +------------------------------- +Constraints and Biased Sampling +------------------------------- + +Source constraints restrict sites to specified cells, materials, or universes, +to fissionable material, or to specified energy and time bounds. Positions +must also lie within the model geometry. For an acceptance region +:math:`\mathcal{A}`, resampling rejected candidates produces the conditional +distribution + +.. math:: + + p_{\mathcal{A}}(X) = + \frac{p(X)\,\mathbf{1}_{\mathcal{A}}(X)} + {\int_{\mathcal{A}}p(Y)\,dY}, + +where :math:`X` denotes the sampled phase-space variables. The base source +wrapper repeats sampling until a candidate passes, or, with the ``kill`` +strategy, returns a failing candidate with weight zero. Independent sources +handle constraints internally: they resample positions for spatial constraints +and energies for energy constraints and the available energy range. Their time +samples are drawn directly. Constraints on an element source in a +:class:`openmc.MeshSource` are applied after selecting that element; constraints +on the mesh source itself are checked by its outer wrapper. Similarly, selecting +among multiple top-level sources precedes each source's constraint handling. +Resampling therefore conditions each selected source separately and retains its +selection probability. + +OpenMC limits the number of consecutive rejections and monitors rejected sites +to detect source definitions that cannot efficiently produce acceptable +particles. The geometry and source domains should be chosen to keep acceptance +reasonable. + +Source biasing draws from a proposal distribution :math:`g` in place of the +desired density :math:`p`, with the correction +:math:`w\leftarrow w\,p(X)/g(X)`. The proposal must cover every region where +the desired distribution contributes. When multiple independent factors are +biased, their corrections multiply. Mesh and point-cloud distributions can +similarly change their discrete selection probabilities and apply the ratio +of the original to the biased normalized probability. The derivation and +available biasing mechanisms are described in :ref:`methods_source_biasing`. diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index bf2ab4ebbb6..387a7602054 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -167,6 +167,9 @@ least 10,000 particles per generation to avoid any bias in the estimate of External Source Distributions ----------------------------- +The sampling algorithms and mathematical derivations for the source classes +are described in :ref:`methods_source_sampling`. + External source distributions can be specified through the :attr:`Settings.source` attribute. If you have a single external source, you can create an instance of any of the subclasses of :class:`openmc.SourceBase` @@ -393,6 +396,9 @@ A surface source file is read back like any other source file, with the Tokamak Plasma Sources ---------------------- +The volume Jacobian, radial distribution, and conditional poloidal sampling +are derived in :ref:`methods_tokamak_source`. + 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 @@ -442,6 +448,9 @@ are born at :math:`t=0`. The toroidal extent can be restricted with Stellarator Plasma Sources -------------------------- +The Fourier volume Jacobian and the radial and angular rejection algorithms +are derived in :ref:`methods_stellarator_source`. + For fully three-dimensional plasmas, the :class:`openmc.StellaratorSource` class samples neutron positions from the flux-surface Fourier representation shared by the `VMEC `_ and `DESC From e6d6e0d02b207a2c5c25aa1c73fa2d8e84653132 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Oct 2026 12:25:01 -0500 Subject: [PATCH 04/15] Clarify term in derivation --- docs/source/methods/source_sampling.rst | 41 +++++++++++++++++++------ 1 file changed, 31 insertions(+), 10 deletions(-) diff --git a/docs/source/methods/source_sampling.rst b/docs/source/methods/source_sampling.rst index fa92a4f13a8..347f71c67b3 100644 --- a/docs/source/methods/source_sampling.rst +++ b/docs/source/methods/source_sampling.rst @@ -477,14 +477,28 @@ Because Fourier evaluation is linear in the coefficients, at fixed angles Z_\rho=\frac{Z_{i+1}-Z_i}{h_i} are constant within the interval. It follows that :math:`\tau` is linear in -:math:`\eta` and :math:`f=\sigma R\tau` is quadratic. Angular integration -preserves this polynomial degree, so the differential volume within interval -:math:`i` is also quadratic: +:math:`\eta` and :math:`f=\sigma R\tau` is quadratic. The differential volume +is the angular integral of this density. Denote its value at +:math:`\rho=\rho_i+h_i\eta` by :math:`V'_i(\eta)`: .. math:: :label: stellarator-source-volume-polynomial - V'_i(\eta)=A_0+A_1\eta+A_2\eta^2. + \begin{aligned} + V'_i(\eta) + &\equiv \left.\frac{dV}{d\rho}\right|_{\rho=\rho_i+h_i\eta} \\ + &= \int_0^{2\pi}\!\int_0^{2\pi} + f_i(\eta,\theta,\zeta)\,d\theta\,d\zeta \\ + &= A_0+A_1\eta+A_2\eta^2. + \end{aligned} + +Angular integration preserves the quadratic dependence on :math:`\eta`, so +the integral and the polynomial are two expressions for the same function. +The coefficients :math:`A_0`, :math:`A_1`, and :math:`A_2` are constant within +each interval. The prime denotes differentiation with respect to :math:`\rho`; +:math:`\eta` specifies where that derivative is evaluated. Since +:math:`d\rho=h_i\,d\eta`, the corresponding volume element is +:math:`dV=h_i V'_i(\eta)\,d\eta`. The derivative generally has different one-sided values at an interior radial surface because the coefficient slopes change there. Each interval uses its @@ -494,8 +508,10 @@ adjacent intervals are kept separate. Computing Differential Volumes ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Initialization evaluates :math:`f` at :math:`\eta=0,1/2,1` on a uniform -angular grid. Its sizes are +To determine the three polynomial coefficients in +:eq:`stellarator-source-volume-polynomial`, initialization computes the angular +integral at just three radial positions, :math:`\eta=0,1/2,1`. At each of +these positions, it evaluates :math:`f_i` on a uniform angular grid with sizes .. math:: @@ -503,7 +519,8 @@ angular grid. Its sizes are N_\zeta=\max(64,8|n|_{\max}). The poloidal grid covers :math:`[0,2\pi)` and the toroidal grid covers one field -period, :math:`[0,2\pi/N_{\mathrm{fp}})`. The full-torus integral is +period, :math:`[0,2\pi/N_{\mathrm{fp}})`. The angular integral defining +:math:`V'_i` is evaluated by the quadrature sum .. math:: @@ -517,8 +534,8 @@ and toroidal degrees at most three times the respective geometry mode limits. The trapezoidal grid therefore integrates this finite Fourier series exactly in exact arithmetic, provided the Jacobian has a consistent sign. -Let :math:`v_0`, :math:`v_m`, and :math:`v_1` be these integrals at -:math:`\eta=0,1/2,1`. Interpolation of the quadratic gives +Let :math:`v_0=V'_i(0)`, :math:`v_m=V'_i(1/2)`, and :math:`v_1=V'_i(1)` be +the three values obtained by quadrature. These determine the coefficients: .. math:: @@ -526,6 +543,9 @@ Let :math:`v_0`, :math:`v_m`, and :math:`v_1` be these integrals at A_1=4v_m-3v_0-v_1, \qquad A_2=2(v_0+v_1-2v_m). +The resulting polynomial supplies :math:`V'_i(\eta)` at any radius within the +interval without further angular integration. + Only two adjacent surfaces' angular grids are held in memory at once. Their positions and poloidal derivatives give the interval's radial derivatives by finite differences. Initialization retains the per-interval integrals and @@ -550,7 +570,8 @@ PDF in the local coordinate: c_3 &= \Delta S_i A_2, \end{aligned} -where :math:`\Delta S_i=S_{i+1}-S_i`. Its integrated mass is +where :math:`\Delta S_i=S_{i+1}-S_i`. Its integrated mass includes the change +of variable :math:`d\rho=h_i\,d\eta`: .. math:: :label: stellarator-source-interval-mass From de3b4788eb184fd00a530687b40ea44fa2d7bdaf Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Oct 2026 12:29:53 -0500 Subject: [PATCH 05/15] Mention discrete Fourier orthogonality --- docs/source/methods/source_sampling.rst | 38 ++++++++++++++++++++++--- 1 file changed, 34 insertions(+), 4 deletions(-) diff --git a/docs/source/methods/source_sampling.rst b/docs/source/methods/source_sampling.rst index 347f71c67b3..a48c6c74d5f 100644 --- a/docs/source/methods/source_sampling.rst +++ b/docs/source/methods/source_sampling.rst @@ -529,10 +529,35 @@ period, :math:`[0,2\pi/N_{\mathrm{fp}})`. The angular integral defining \sum_{j=0}^{N_\theta-1}\sum_{k=0}^{N_\zeta-1} f_i(\eta,\theta_j,\zeta_k). -The weight includes all field periods. Products in :math:`R\tau` have poloidal -and toroidal degrees at most three times the respective geometry mode limits. -The trapezoidal grid therefore integrates this finite Fourier series exactly -in exact arithmetic, provided the Jacobian has a consistent sign. +The weight includes all field periods. Exactness of this quadrature follows +from discrete Fourier orthogonality, or the roots-of-unity identity: + +.. math:: + + \frac{1}{N}\sum_{j=0}^{N-1} e^{2\pi\mathrm{i}\ell j/N} + = \begin{cases} + 1, & \ell \text{ is a multiple of } N, \\ + 0, & \text{otherwise}. + \end{cases} + +This identity follows from summing a finite geometric series. For a +trigonometric polynomial with modes :math:`-K,\ldots,K`, choosing +:math:`N>K` makes every nonconstant mode sum to zero. The periodic trapezoidal +rule then returns exactly the integral of the constant mode. This result is +stated in Corollary 3.3 of [TrefethenWeideman]_. + +Products in :math:`R\tau` have poloidal and toroidal degrees at most three +times the respective geometry mode limits. Applying the identity independently +in both angles gives sufficient conditions for exact integration: + +.. math:: + + N_\theta>3m_{\max}, \qquad N_\zeta>3|n|_{\max}. + +The grid sizes above satisfy these conditions. The quadrature is thus exact +in exact arithmetic, provided the Jacobian has a consistent sign. Recovering +only the integral requires fewer points than reconstructing the full Fourier +series under the Nyquist sampling condition. Let :math:`v_0=V'_i(0)`, :math:`v_m=V'_i(1/2)`, and :math:`v_1=V'_i(1)` be the three values obtained by quadrature. These determine the coefficients: @@ -739,3 +764,8 @@ biased, their corrections multiply. Mesh and point-cloud distributions can similarly change their discrete selection probabilities and apply the ratio of the original to the biased normalized probability. The derivation and available biasing mechanisms are described in :ref:`methods_source_biasing`. + +.. [TrefethenWeideman] Lloyd N. Trefethen and J. A. C. Weideman, + `The Exponentially Convergent Trapezoidal Rule + `_, + *SIAM Review*, **56** (3), 385--458 (2014), Corollary 3.3, p. 398. From dea579304cf7f82cba0c9851c906bc0669d4582a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 21 Oct 2026 09:46:13 -0500 Subject: [PATCH 06/15] Update to reflect new angular sampling strategy --- docs/source/methods/source_sampling.rst | 198 ++++++++++++++++++++---- docs/source/usersguide/settings.rst | 9 +- 2 files changed, 177 insertions(+), 30 deletions(-) diff --git a/docs/source/methods/source_sampling.rst b/docs/source/methods/source_sampling.rst index a48c6c74d5f..ac427e617e6 100644 --- a/docs/source/methods/source_sampling.rst +++ b/docs/source/methods/source_sampling.rst @@ -649,8 +649,10 @@ proportion to its acceptance rate. After acceptance, :math:`\rho=\rho_i+h_i\eta` has the required marginal :math:`p_\rho\propto S V'`. -Conditional Angular Rejection -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +.. _methods_stellarator_angular_sampling: + +Conditional Angular Sampling +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ With :math:`\rho` fixed, the target PDF is @@ -659,36 +661,176 @@ With :math:`\rho` fixed, the target PDF is p(\theta,\zeta\mid\rho) = \frac{f_i(\eta,\theta,\zeta)}{V'_i(\eta)}. -At each angular grid point, the three radial evaluations determine the -quadratic :math:`f_i(\eta)` throughout the interval. Its maximum is obtained -from the endpoints and, for a concave quadratic, an interior vertex. The current -angular envelope is +The angular proposal uses a 32 by 32 partition of one field period, with a +separate bound for each cell that varies with radius. To describe the cell +geometry, introduce the field-period angle +:math:`\varphi=N_{\mathrm{fp}}\zeta` and write + +.. math:: + + F_i(\eta,\theta,\varphi) + = f_i(\eta,\theta,\varphi/N_{\mathrm{fp}}). + +Both :math:`\theta` and :math:`\varphi` range over :math:`[0,2\pi)`, so each +cell has angular half-width :math:`d=\pi/32` in both coordinates. This cell grid +is separate from the finer quadrature grid used during initialization. The +quadrature grid resolves the density's Fourier coefficients; the cell grid +defines regions for proposing samples. Integrating :math:`F_i` over the +rescaled square gives :math:`V'_i(\eta)`: the factor +:math:`d\varphi=N_{\mathrm{fp}}\,d\zeta` accounts for all field periods. + +Within a fixed radial interval, express the quadratic dependence on +:math:`\eta` in the Bernstein basis: + +.. math:: + :label: stellarator-source-angular-bernstein + + \begin{aligned} + F_i(\eta,\theta,\varphi) + &= \sum_{j=0}^2 b_j(\eta)G_j(\theta,\varphi), \\ + b_0(\eta)&=(1-\eta)^2, \qquad + b_1(\eta)=2\eta(1-\eta), \qquad b_2(\eta)=\eta^2. + \end{aligned} + +Suppressing the interval index on the three angular control functions, their +values follow from the radial evaluations already available during setup: + +.. math:: + + \begin{aligned} + G_0 &= F_i(0,\theta,\varphi), \\ + G_1 &= 2F_i(1/2,\theta,\varphi) + -\tfrac12[F_i(0,\theta,\varphi)+F_i(1,\theta,\varphi)], \\ + G_2 &= F_i(1,\theta,\varphi). + \end{aligned} + +This is an exact change of polynomial basis. The useful property is that +:math:`b_j\ge0` and :math:`\sum_j b_j=1` on the interval. If a nonnegative +constant :math:`U_{jc}` bounds :math:`G_j` throughout angular cell :math:`c`, +then + +.. math:: + :label: stellarator-source-cell-bound + + B_c(\eta)=\sum_{j=0}^2 b_j(\eta)U_{jc} + \;\ge\; F_i(\eta,\theta,\varphi) + \quad\text{within cell }c. + +Thus three bounds computed once per cell produce a bound at every radius in +the interval. The control functions :math:`G_j` can have negative values even +when the density is nonnegative; they are bounded functions, not probability +densities. Nonnegativity of the Bernstein basis weights is what preserves the +upper-bound inequality. + +Bounding the Density on Angular Cells +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Each :math:`G_j` is a finite Fourier series in the two angular coordinates: + +.. math:: + + G_j(\theta,\varphi) + = \sum_{m=-M}^{M}\sum_{n=-L}^{L} + \widehat G_{j,mn}e^{\mathrm{i}(m\theta+n\varphi)}, + \qquad M=3m_{\max},\quad L=3|n|_{\max}. + +Initialization recovers these coefficients using separable discrete Fourier +sums on the quadrature grid. Resolving all coefficients requires the stronger +conditions :math:`N_\theta>2M` and :math:`N_\zeta>2L`, which the grid sizes +above also satisfy. The implementation evaluates the sums directly and needs +no external FFT library. + +Let :math:`\boldsymbol{a}_c=(\theta_c,\varphi_c)` be a cell center. Any point +in the cell has displacement :math:`\boldsymbol{v}` with +:math:`|v_\theta|,|v_\varphi|\le d`. Taylor's theorem gives the upper bound + +.. math:: + :label: stellarator-source-cell-taylor + + \begin{aligned} + G_j(\boldsymbol{a}_c+\boldsymbol{v}) + &\le G_j(\boldsymbol{a}_c) + +d\left(|\partial_\theta G_j(\boldsymbol{a}_c)| + +|\partial_\varphi G_j(\boldsymbol{a}_c)|\right) \\ + &\qquad +\frac{d^2}{2}\sum_{m,n} + |\widehat G_{j,mn}|(|m|+|n|)^2. + \end{aligned} + +To obtain the remainder term, consider the second derivative along the line +from the center to the proposed point. Each Fourier term contributes a factor +:math:`(mv_\theta+nv_\varphi)^2`, bounded by +:math:`d^2(|m|+|n|)^2`. Summing absolute Fourier amplitudes bounds that +derivative everywhere along the line. The center value and its first +derivatives retain the cancellations among Fourier modes, allowing a tighter +bound than separately bounding all factors of the Jacobian. + +OpenMC adds a roundoff allowance of +:math:`10^{-10}\sum_{m,n}|\widehat G_{j,mn}|` to the right-hand side of +:eq:`stellarator-source-cell-taylor`, then takes its maximum with zero to form +:math:`U_{jc}`. Center values and derivatives are evaluated using all recovered +Fourier coefficients, including frequencies above the cell grid's Nyquist +limit. They are folded onto the cell-center grid with their phase and +derivative factors preserved; the cell grid does not truncate the density. +Bounds for both Jacobian orientations share these derivative bounds, with the +center value changing sign. Initialization retains the bounds for the +orientation determined by the total signed volume. + +Selecting and Sampling Angular Cells +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +For each control function, define the sum of its cell bounds + +.. math:: + + T_j=\sum_c U_{jc}. + +Initialization builds three alias tables, one for each distribution +:math:`U_{jc}/T_j`. At the sampled radius, OpenMC selects a control index +:math:`j` with probability + +.. math:: + + P(j\mid\eta)=\frac{b_j(\eta)T_j}{\sum_k b_k(\eta)T_k}, + +then selects a cell from that control's alias table. Summing over control +indices shows that the resulting cell probability is + +.. math:: + + P(c\mid\eta) + =\sum_{j:T_j>0} P(j\mid\eta)\frac{U_{jc}}{T_j} + =\frac{B_c(\eta)}{\sum_{c'}B_{c'}(\eta)}. + +Components with :math:`T_j=0` have zero selection probability and require no +alias table. This mixture samples the radius-dependent bounds using three +fixed tables, avoiding construction of a new table for each sampled radius. + +The sampler proposes a point uniformly inside the chosen cell and accepts it +when .. math:: - B_i^{\mathrm{ang}} = 1.05\max_{j,k}\max_{0\le\eta\le1} - f_i(\eta,\theta_j,\zeta_k). + \xi B_c(\eta) < F_i(\eta,\theta,\varphi). -The radial maximization here is exact in exact arithmetic. The angular maximum -is estimated on a grid, with a 5% heuristic allowance for peaks between nodes. -Exact Fourier quadrature of an integral does not establish an upper bound on -the function. This envelope therefore does not provide a rigorous bound for -every supported coefficient table; correctness of angular rejection requires -it to bound the continuous angular density. +All cells have the same area, so the proposal density within a cell is +proportional to :math:`B_c`. Multiplying it by the acceptance probability +:math:`F_i/B_c` gives a density proportional to :math:`F_i`, as required. +On rejection, both the control index and the cell are redrawn while retaining +the radial interval and :math:`\rho`. -The sampler proposes :math:`\theta=2\pi\xi_1` and -:math:`\zeta=2\pi\xi_2` over the full torus, and accepts when +After acceptance, an independently and uniformly chosen field-period index +:math:`\ell\in\{0,\ldots,N_{\mathrm{fp}}-1\}` gives the full-torus angle .. math:: - \xi_3 B_i^{\mathrm{ang}} < f_i(\eta,\theta,\zeta). + \zeta=\frac{\varphi+2\pi\ell}{N_{\mathrm{fp}}}. -Rejected angles are redrawn while retaining both the radial interval and -:math:`\rho`. If the envelope is valid, this produces the angular conditional -PDF above. An underestimated envelope clips the peaks because acceptance -cannot exceed one; the current sampling loop has no explicit check for that -condition. The angular rejection limit detects excessive rejection, which is a -different failure condition. +The sampled :math:`R` and :math:`Z` apply in every field period. The cell +partition affects rejection efficiency while the accepted points follow the +continuous Fourier density. The Taylor inequality supplies bounds in exact +arithmetic; a runtime check terminates sampling if a density exceeds its cell +bound. Nonfinite or negative densities and excessive rejection also produce +errors. Sampling Sequence and Approximation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -698,8 +840,9 @@ The steps in the stellarator source's sampling method are: 1. Select a radial interval using its integrated mass :math:`M_i`. 2. Sample the local radial coordinate by cubic-PDF rejection within that interval, and calculate :math:`\rho`. -3. Sample the poloidal and toroidal angles by conditional rejection at that - fixed radius. +3. At that fixed radius, select a control index and an angular cell, propose + angles uniformly inside the cell, and reject using the local bound. Repeat + angular proposals until accepted, then choose a field period uniformly. 4. Map the accepted position to :math:`(R\cos\zeta,R\sin\zeta,Z)`. 5. Sample an isotropic flight direction. 6. Sample energy from the common spectrum or the mixture of spectra bracketing @@ -709,8 +852,9 @@ Radial integration and rejection sampling are consistent with the piecewise linear coefficient and emission tables, up to floating-point effects. Those tables approximate the underlying equilibrium and emission profile; increasing the number of radial surfaces can reduce interpolation error. The finite -Fourier representation also limits the geometric detail retained. These input -approximations are distinct from the angular-envelope condition above. +Fourier representation also limits the geometric detail retained. The +Bernstein representation and angular-cell partition introduce no additional +approximation to this density; the angular bounds control rejection efficiency. The :meth:`~openmc.StellaratorSource.from_vmec` and :meth:`~openmc.StellaratorSource.from_desc` constructors prepare the coefficient diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 387a7602054..16e45cf579b 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -485,12 +485,15 @@ radial proposals are resampled within the same interval. The two angles are sampled from the conditional distribution :math:`p(\theta,\zeta|\rho) \propto R\,|\tau|` (with :math:`\tau` the poloidal-plane Jacobian) by rejection. -Angular proposals use a 32 by 32 partition of one field period. Each cell has -a bound that varies quadratically across the radial interval. OpenMC selects +Angular proposals use a 32 by 32 partition of one field period in the two +spatial angular coordinates. Each cell has a bound that varies quadratically +across the radial interval. OpenMC selects cells in proportion to these bounds, samples uniformly inside the chosen cell, and accepts with probability equal to the density divided by the cell bound. Rejected angular proposals retain the radius and select a new cell. A uniformly -chosen field period maps the accepted point to the full torus. +chosen field period maps the accepted point to the full torus. The construction +of the bounds and the Bernstein mixture used to select cells are derived in +:ref:`methods_stellarator_angular_sampling`. The bounds are constructed during C++ initialization from the density's Fourier coefficients, using a Taylor bound on each cell. This requires no additional From 430244a5517d70ad57697329bb5fadddd38bd6f7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 21 Oct 2026 22:29:53 -0500 Subject: [PATCH 07/15] Remove incorrect note --- docs/source/methods/source_sampling.rst | 7 ------- 1 file changed, 7 deletions(-) diff --git a/docs/source/methods/source_sampling.rst b/docs/source/methods/source_sampling.rst index ac427e617e6..48421c55e52 100644 --- a/docs/source/methods/source_sampling.rst +++ b/docs/source/methods/source_sampling.rst @@ -396,13 +396,6 @@ The position is then mapped to :math:`(R\cos\phi,R\sin\phi,Z)`. This construction samples the spatial coordinates without a rejection loop, with discretization error from the radial and poloidal tabulations. -.. note:: - - The continuous integrals above give the required mixture weights. The - current implementation stores different values for :math:`I_2`, - :math:`I_4`, and :math:`I_5`; that mismatch introduces an additional error - in poloidal sampling beyond the tabulation error. - .. _methods_stellarator_source: Stellarator Sources From 49e832331b9180556bc4f9f60e4e193445bb1871 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 24 Oct 2026 08:18:15 -0500 Subject: [PATCH 08/15] Better derivation of angular density bounds --- docs/source/methods/source_sampling.rst | 104 +++++++++++++++++++++--- 1 file changed, 93 insertions(+), 11 deletions(-) diff --git a/docs/source/methods/source_sampling.rst b/docs/source/methods/source_sampling.rst index 48421c55e52..52a901d3579 100644 --- a/docs/source/methods/source_sampling.rst +++ b/docs/source/methods/source_sampling.rst @@ -734,8 +734,90 @@ above also satisfy. The implementation evaluates the sums directly and needs no external FFT library. Let :math:`\boldsymbol{a}_c=(\theta_c,\varphi_c)` be a cell center. Any point -in the cell has displacement :math:`\boldsymbol{v}` with -:math:`|v_\theta|,|v_\varphi|\le d`. Taylor's theorem gives the upper bound +in the cell can be written as :math:`\boldsymbol{a}_c+\boldsymbol{v}`, where +:math:`|v_\theta|,|v_\varphi|\le d`. The goal is to bound the largest possible +increase from :math:`G_j(\boldsymbol{a}_c)` over all such displacements. Taylor's +theorem separates this increase into a contribution from the slope at the +center and a remainder accounting for curvature. + +For a fixed displacement, reduce the problem to one dimension by following +the straight line from the center to the point: + +.. math:: + + g(\lambda)=G_j(\boldsymbol{a}_c+\lambda\boldsymbol{v}), + \qquad 0\le\lambda\le1. + +Here :math:`\lambda` measures progress along the line; the radial coordinate +:math:`\eta` remains fixed during angular sampling. Taylor's theorem with the +Lagrange remainder gives + +.. math:: + :label: stellarator-source-line-taylor + + g(1)=g(0)+g'(0)+\frac12 g''(\lambda_*), + \qquad \lambda_*\in(0,1). + +This is an exact identity. The second derivative is evaluated at an unknown +intermediate point, and its term accounts for the entire departure from the +linear approximation. Bounding this derivative along the line therefore bounds +the full remainder, including effects that a longer Taylor expansion would +express through higher derivatives. No additional terms need to be neglected. + +The center value and slope follow from the chain rule: + +.. math:: + + \begin{aligned} + g(0)&=G_j(\boldsymbol{a}_c), \\ + g'(0)&=v_\theta\partial_\theta G_j(\boldsymbol{a}_c) + +v_\varphi\partial_\varphi G_j(\boldsymbol{a}_c). + \end{aligned} + +Using the maximum displacement in each coordinate bounds the slope term by + +.. math:: + + g'(0)\le d\left(|\partial_\theta G_j(\boldsymbol{a}_c)| + +|\partial_\varphi G_j(\boldsymbol{a}_c)|\right). + +To bound the remainder, differentiate the Fourier series twice along the line: + +.. math:: + + \begin{aligned} + g''(\lambda) + &= -\sum_{m,n}\widehat G_{j,mn}(mv_\theta+nv_\varphi)^2 \\ + &\qquad\times e^{\mathrm{i}[m(\theta_c+\lambda v_\theta) + +n(\varphi_c+\lambda v_\varphi)]}. + \end{aligned} + +Each exponential has magnitude one. The triangle inequality then gives a +bound valid at every point along the line: + +.. math:: + + |g''(\lambda)|\le\sum_{m,n}|\widehat G_{j,mn}| + |mv_\theta+nv_\varphi|^2. + +Because the point lies within the cell, + +.. math:: + + |mv_\theta+nv_\varphi| + \le |m|\,|v_\theta|+|n|\,|v_\varphi| + \le d(|m|+|n|). + +Consequently, + +.. math:: + + |g''(\lambda)|\le d^2 C_j, \qquad + C_j=\sum_{m,n}|\widehat G_{j,mn}|(|m|+|n|)^2. + +The quantity :math:`C_j` depends only on the Fourier coefficients and bounds +curvature for every cell and displacement. Substituting this result and the +slope bound into :eq:`stellarator-source-line-taylor` gives .. math:: :label: stellarator-source-cell-taylor @@ -745,17 +827,17 @@ in the cell has displacement :math:`\boldsymbol{v}` with &\le G_j(\boldsymbol{a}_c) +d\left(|\partial_\theta G_j(\boldsymbol{a}_c)| +|\partial_\varphi G_j(\boldsymbol{a}_c)|\right) \\ - &\qquad +\frac{d^2}{2}\sum_{m,n} - |\widehat G_{j,mn}|(|m|+|n|)^2. + &\qquad +\frac{d^2}{2}C_j. \end{aligned} -To obtain the remainder term, consider the second derivative along the line -from the center to the proposed point. Each Fourier term contributes a factor -:math:`(mv_\theta+nv_\varphi)^2`, bounded by -:math:`d^2(|m|+|n|)^2`. Summing absolute Fourier amplitudes bounds that -derivative everywhere along the line. The center value and its first -derivatives retain the cancellations among Fourier modes, allowing a tighter -bound than separately bounding all factors of the Jacobian. +Every quantity on the right-hand side is available during initialization and +is independent of the unknown displacement. The inequality therefore bounds +the function throughout the cell, including peaks between evaluation points. +The curvature allowance remains present even when the center gradient is zero. +The center value and its first derivatives retain the cancellations among +Fourier modes, allowing a tighter bound than separately bounding all factors of +the Jacobian. Taking absolute values in the slope and remainder terms makes +the bound conservative; it need not equal a value attained inside the cell. OpenMC adds a roundoff allowance of :math:`10^{-10}\sum_{m,n}|\widehat G_{j,mn}|` to the right-hand side of From 1f8c41c9b1bafdc6c4418e712fab7f69177590e6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 24 Oct 2026 08:26:43 -0500 Subject: [PATCH 09/15] Simplify last paragraph --- docs/source/methods/source_sampling.rst | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/docs/source/methods/source_sampling.rst b/docs/source/methods/source_sampling.rst index 52a901d3579..db794f0b696 100644 --- a/docs/source/methods/source_sampling.rst +++ b/docs/source/methods/source_sampling.rst @@ -839,12 +839,10 @@ Fourier modes, allowing a tighter bound than separately bounding all factors of the Jacobian. Taking absolute values in the slope and remainder terms makes the bound conservative; it need not equal a value attained inside the cell. -OpenMC adds a roundoff allowance of -:math:`10^{-10}\sum_{m,n}|\widehat G_{j,mn}|` to the right-hand side of -:eq:`stellarator-source-cell-taylor`, then takes its maximum with zero to form -:math:`U_{jc}`. Center values and derivatives are evaluated using all recovered -Fourier coefficients, including frequencies above the cell grid's Nyquist -limit. They are folded onto the cell-center grid with their phase and +Taking the maximum of the right-hand side of +:eq:`stellarator-source-cell-taylor` gives a nonnegative bound :math:`U_{jc}`. +Center values and derivatives are evaluated using all recovered Fourier +coefficients. They are folded onto the cell-center grid with their phase and derivative factors preserved; the cell grid does not truncate the density. Bounds for both Jacobian orientations share these derivative bounds, with the center value changing sign. Initialization retains the bounds for the From 242f683d85d18f21260e9e9055e7856de0f3bbad Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 24 Oct 2026 08:43:21 -0500 Subject: [PATCH 10/15] More improvements --- docs/source/methods/source_sampling.rst | 35 ++++++++++++++----------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/docs/source/methods/source_sampling.rst b/docs/source/methods/source_sampling.rst index db794f0b696..256f3d4acf9 100644 --- a/docs/source/methods/source_sampling.rst +++ b/docs/source/methods/source_sampling.rst @@ -440,13 +440,17 @@ and :math:`Z_\zeta` drop out of the determinant because the cylindrical azimuth is exactly :math:`\zeta`. Thus the joint coordinate PDF is :math:`p(\rho,\theta,\zeta)\propto S(\rho)R|\tau|`. -The map is assumed to have positive :math:`R`, nested surfaces, and a -consistently oriented Jacobian away from the magnetic axis. Either orientation -is allowed. Initialization determines a global sign :math:`\sigma` from a -signed volume estimate and evaluates :math:`f=\sigma R\tau`, so that -:math:`f=R|\tau|` for a valid map. Grid checks reject nonpositive radii or -significant changes in Jacobian sign; these checks are not a global geometric -certification. +The geometry must have positive :math:`R` and nested flux surfaces. The sign +of :math:`\tau` depends on the convention for measuring the poloidal angle: +traversing a cross section clockwise or counterclockwise describes the same +surface but gives opposite signs for :math:`\tau`. Physical volume uses its +absolute value, so either convention is allowed. For a valid map, the sign is +consistent throughout the plasma away from the magnetic axis, where the +Jacobian vanishes. OpenMC determines this sign during initialization and sets +:math:`\sigma=+1` for positive :math:`\tau` or :math:`\sigma=-1` for negative +:math:`\tau`. It can then evaluate the positive sampling density as +:math:`f=\sigma R\tau=R|\tau|`. Initialization checks for nonpositive radii +and significant changes in the sign of :math:`\tau` on a setup grid. Interpolation Within a Radial Interval ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -839,14 +843,15 @@ Fourier modes, allowing a tighter bound than separately bounding all factors of the Jacobian. Taking absolute values in the slope and remainder terms makes the bound conservative; it need not equal a value attained inside the cell. -Taking the maximum of the right-hand side of -:eq:`stellarator-source-cell-taylor` gives a nonnegative bound :math:`U_{jc}`. -Center values and derivatives are evaluated using all recovered Fourier -coefficients. They are folded onto the cell-center grid with their phase and -derivative factors preserved; the cell grid does not truncate the density. -Bounds for both Jacobian orientations share these derivative bounds, with the -center value changing sign. Initialization retains the bounds for the -orientation determined by the total signed volume. +For each control function :math:`G_j` and cell :math:`c`, OpenMC evaluates the +right-hand side of :eq:`stellarator-source-cell-taylor` using the function's +value and slopes at the cell center and the curvature bound :math:`C_j`. +All Fourier coefficients contribute to these calculations, so the bound +accounts for variation inside the cell as well as at its center. If the result +is negative, OpenMC uses zero instead, giving a nonnegative bound +:math:`U_{jc}` suitable for constructing the cell-selection probabilities. +These bounds are stored during initialization and combined at the sampled +radius as described below. Selecting and Sampling Angular Cells ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ From ae4ac971c88bf822956c18ffed5dae247fc53687 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 24 Oct 2026 08:50:34 -0500 Subject: [PATCH 11/15] Be clear about nonnegativity --- docs/source/methods/source_sampling.rst | 31 ++++++++++++++++++------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/docs/source/methods/source_sampling.rst b/docs/source/methods/source_sampling.rst index 256f3d4acf9..6498ef18007 100644 --- a/docs/source/methods/source_sampling.rst +++ b/docs/source/methods/source_sampling.rst @@ -714,10 +714,25 @@ then \quad\text{within cell }c. Thus three bounds computed once per cell produce a bound at every radius in -the interval. The control functions :math:`G_j` can have negative values even -when the density is nonnegative; they are bounded functions, not probability -densities. Nonnegativity of the Bernstein basis weights is what preserves the -upper-bound inequality. +the interval. For a valid geometry, the control functions :math:`G_j` are +also nonnegative. This follows from the factorization +:math:`F_i=R\widetilde\tau`, where :math:`\widetilde\tau=\sigma\tau` is the +nonnegative poloidal Jacobian. Both factors are linear in :math:`\eta`. +Writing their endpoint values as :math:`R_0,R_1` and +:math:`\widetilde\tau_0,\widetilde\tau_1` gives + +.. math:: + + \begin{aligned} + G_0 &= R_0\widetilde\tau_0, \\ + G_1 &= \tfrac12(R_0\widetilde\tau_1+R_1\widetilde\tau_0), \\ + G_2 &= R_1\widetilde\tau_1. + \end{aligned} + +Each expression is nonnegative because the endpoint radii are positive and +the endpoint values of :math:`\widetilde\tau` are nonnegative. The +:math:`G_j` are unnormalized angular functions; the cell-selection +probabilities are constructed from their bounds below. Bounding the Density on Angular Cells ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -847,11 +862,9 @@ For each control function :math:`G_j` and cell :math:`c`, OpenMC evaluates the right-hand side of :eq:`stellarator-source-cell-taylor` using the function's value and slopes at the cell center and the curvature bound :math:`C_j`. All Fourier coefficients contribute to these calculations, so the bound -accounts for variation inside the cell as well as at its center. If the result -is negative, OpenMC uses zero instead, giving a nonnegative bound -:math:`U_{jc}` suitable for constructing the cell-selection probabilities. -These bounds are stored during initialization and combined at the sampled -radius as described below. +accounts for variation inside the cell as well as at its center. The resulting +bounds :math:`U_{jc}` are stored during initialization and combined at the +sampled radius as described below. Selecting and Sampling Angular Cells ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ From 018d5746221dde878eb43acdb112ef63c2ab8722 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 24 Oct 2026 10:19:55 -0500 Subject: [PATCH 12/15] Consistency of notation in source sampling docs --- docs/source/methods/source_sampling.rst | 192 +++++++++++++++--------- 1 file changed, 123 insertions(+), 69 deletions(-) diff --git a/docs/source/methods/source_sampling.rst b/docs/source/methods/source_sampling.rst index 6498ef18007..a1c0ceb0904 100644 --- a/docs/source/methods/source_sampling.rst +++ b/docs/source/methods/source_sampling.rst @@ -16,8 +16,11 @@ This chapter describes the algorithms represented by the subclasses of :class:`openmc.SourceBase`. Input examples are given in :ref:`usersguide_source`. Unless otherwise stated, :math:`\xi` denotes an independent uniform random number on :math:`[0,1)`, generated as described in -:ref:`methods_random_numbers`. Probability densities below are normalized with -respect to the coordinates explicitly indicated. +:ref:`methods_random_numbers`. Normalized probability densities are denoted by +:math:`p` or :math:`f`, cumulative distribution functions by :math:`F`, and +discrete selection probabilities by :math:`P`. A tilde on :math:`p` indicates +an unnormalized density. Geometric volume factors are denoted by :math:`J`. +Densities are defined with respect to the coordinates explicitly indicated. .. _methods_source_selection: @@ -75,8 +78,9 @@ be continuous, discrete, or a mixture of distributions. Sampling One-Dimensional Distributions -------------------------------------- -For a normalized density :math:`f(x)`, inverse-transform sampling draws -:math:`x=F^{-1}(\xi)`, where +For a scalar random variable :math:`X` with cumulative distribution function +:math:`F(x)=\Pr(X\le x)`, inverse-transform sampling draws +:math:`x=F^{-1}(\xi)`. If the variable has a normalized density :math:`f`, then .. math:: @@ -86,7 +90,7 @@ OpenMC provides analytic distributions as well as tabulated distributions through :mod:`openmc.stats`. A tabulated distribution first selects an interval using its cumulative probability, then inverts the CDF within that interval according to the specified interpolation. For example, a linear-linear PDF on -:math:`[x_i,x_{i+1}]` has +:math:`[x_i,x_{i+1}]`, with :math:`f_i=f(x_i)`, has .. math:: @@ -163,13 +167,15 @@ Sources Defined on Meshes ------------------------- A mesh distribution first selects an element, then draws a position uniformly -within that element. If :math:`a_j` is the integrated strength assigned to -element :math:`j`, with volume :math:`V_j`, the spatial PDF is +within that element. Denote the region occupied by element :math:`j` by +:math:`\mathcal D_j` and its volume by +:math:`V_j=\operatorname{vol}(\mathcal D_j)`. If :math:`a_j` is +the element's integrated strength, the spatial PDF is .. math:: p_r(\mathbf{r}) = \frac{a_j}{\sum_k a_k}\frac{1}{V_j}, - \qquad \mathbf{r}\in V_j. + \qquad \mathbf{r}\in\mathcal D_j. An :class:`openmc.IndependentSource` with a :class:`~openmc.stats.MeshSpatial` distribution uses the same direction, energy, @@ -249,8 +255,10 @@ Radially Varying Spectra ------------------------ Both plasma sources accept either one energy distribution or an energy -distribution at every radial grid point. In an interval with local coordinate -:math:`\eta\in[0,1]`, the conditional energy distribution is the mixture +distribution at every radial grid point. Between grid points :math:`u_i` and +:math:`u_{i+1}`, define the local coordinate +:math:`\eta=(u-u_i)/(u_{i+1}-u_i)\in[0,1]`. The conditional energy distribution +is the mixture .. math:: @@ -292,10 +300,11 @@ Differentiating the geometry gives &= \kappa r\,[G(\alpha)-2Dq\cos\alpha], \\ G(\alpha) &= \cos(\delta\sin\alpha) - +\delta\sin\alpha\cos\alpha\sin\psi. + +\delta\sin\alpha\cos\alpha\sin\psi \end{aligned} -The coordinate volume element is :math:`R` times this poloidal Jacobian. +where subscripts :math:`r` and :math:`\alpha` denote partial derivatives. The +coordinate volume element is :math:`R` times this poloidal Jacobian. Up to constants independent of :math:`q` and :math:`\alpha`, the joint density is therefore @@ -408,7 +417,10 @@ The stellarator source uses flux coordinates :math:`(\rho,\theta,\zeta)`, where :math:`\rho=\sqrt{s}` for normalized toroidal flux :math:`s`, :math:`\theta` is a poloidal angle, and :math:`\zeta` equals the cylindrical azimuth. With :math:`N_{\mathrm{fp}}` field periods, define -:math:`\chi_k=m_k\theta-n_kN_{\mathrm{fp}}\zeta`. The surfaces are +:math:`\chi_k=m_k\theta-n_kN_{\mathrm{fp}}\zeta`, where :math:`m_k\ge0` +and :math:`n_k` are the integer poloidal and toroidal mode numbers. Their +maximum magnitudes are :math:`m_{\max}=\max_k m_k` and +:math:`|n|_{\max}=\max_k|n_k|`. The surfaces are .. math:: :label: stellarator-source-fourier @@ -448,8 +460,8 @@ absolute value, so either convention is allowed. For a valid map, the sign is consistent throughout the plasma away from the magnetic axis, where the Jacobian vanishes. OpenMC determines this sign during initialization and sets :math:`\sigma=+1` for positive :math:`\tau` or :math:`\sigma=-1` for negative -:math:`\tau`. It can then evaluate the positive sampling density as -:math:`f=\sigma R\tau=R|\tau|`. Initialization checks for nonpositive radii +:math:`\tau`. It can then evaluate the nonnegative volume factor as +:math:`J=\sigma R\tau=R|\tau|`. Initialization checks for nonpositive radii and significant changes in the sign of :math:`\tau` on a setup grid. Interpolation Within a Radial Interval @@ -457,13 +469,20 @@ Interpolation Within a Radial Interval On interval :math:`i`, define its width :math:`h_i=\rho_{i+1}-\rho_i` and local coordinate :math:`\eta=(\rho-\rho_i)/h_i`. Each Fourier coefficient :math:`C` -and the emission profile are interpolated linearly: +and the emission profile are interpolated linearly, using endpoint values +:math:`C_i=C(\rho_i)` and :math:`S_i=S(\rho_i)`: .. math:: C(\eta)=(1-\eta)C_i+\eta C_{i+1}, \qquad S(\eta)=S_i+\eta(S_{i+1}-S_i). +Functions written with argument :math:`\eta` are evaluated at +:math:`\rho=\rho_i+h_i\eta`. Throughout the following derivation, a single +interval :math:`i` is considered at a time. Polynomial coefficients, control +functions, and cell bounds depend on that interval; the index :math:`i` is +suppressed where it is not written explicitly. + Because Fourier evaluation is linear in the coefficients, at fixed angles :math:`R`, :math:`Z`, :math:`R_\theta`, and :math:`Z_\theta` are linear in :math:`\eta`, while @@ -474,9 +493,15 @@ Because Fourier evaluation is linear in the coefficients, at fixed angles Z_\rho=\frac{Z_{i+1}-Z_i}{h_i} are constant within the interval. It follows that :math:`\tau` is linear in -:math:`\eta` and :math:`f=\sigma R\tau` is quadratic. The differential volume -is the angular integral of this density. Denote its value at -:math:`\rho=\rho_i+h_i\eta` by :math:`V'_i(\eta)`: +:math:`\eta` and :math:`J=\sigma R\tau` is quadratic. Define the local volume +factor explicitly as + +.. math:: + + J_i(\eta,\theta,\zeta)=J(\rho_i+h_i\eta,\theta,\zeta). + +The differential volume is the angular integral of this factor. Denote its +value at :math:`\rho=\rho_i+h_i\eta` by :math:`V'_i(\eta)`: .. math:: :label: stellarator-source-volume-polynomial @@ -485,7 +510,7 @@ is the angular integral of this density. Denote its value at V'_i(\eta) &\equiv \left.\frac{dV}{d\rho}\right|_{\rho=\rho_i+h_i\eta} \\ &= \int_0^{2\pi}\!\int_0^{2\pi} - f_i(\eta,\theta,\zeta)\,d\theta\,d\zeta \\ + J_i(\eta,\theta,\zeta)\,d\theta\,d\zeta \\ &= A_0+A_1\eta+A_2\eta^2. \end{aligned} @@ -508,7 +533,7 @@ Computing Differential Volumes To determine the three polynomial coefficients in :eq:`stellarator-source-volume-polynomial`, initialization computes the angular integral at just three radial positions, :math:`\eta=0,1/2,1`. At each of -these positions, it evaluates :math:`f_i` on a uniform angular grid with sizes +these positions, it evaluates :math:`J_i` on a uniform angular grid with sizes .. math:: @@ -524,7 +549,7 @@ period, :math:`[0,2\pi/N_{\mathrm{fp}})`. The angular integral defining V'_i(\eta) = \frac{4\pi^2}{N_\theta N_\zeta} \sum_{j=0}^{N_\theta-1}\sum_{k=0}^{N_\zeta-1} - f_i(\eta,\theta_j,\zeta_k). + J_i(\eta,\theta_j,\zeta_k). The weight includes all field periods. Exactness of this quadrature follows from discrete Fourier orthogonality, or the roots-of-unity identity: @@ -578,13 +603,13 @@ Selecting and Sampling a Radial Interval Multiplying the linear emission profile by :eq:`stellarator-source-volume-polynomial` yields an unnormalized cubic radial -PDF in the local coordinate: +density with respect to :math:`d\rho`, expressed using the local coordinate: .. math:: :label: stellarator-source-radial-polynomial \begin{aligned} - P_i(\eta) &= S(\eta)V'_i(\eta) + \widetilde p_i(\eta) &= S(\eta)V'_i(\eta) = c_0+c_1\eta+c_2\eta^2+c_3\eta^3, \\ c_0 &= S_i A_0, \\ c_1 &= S_i A_1+\Delta S_i A_0, \\ @@ -599,7 +624,7 @@ of variable :math:`d\rho=h_i\,d\eta`: :label: stellarator-source-interval-mass \begin{aligned} - M_i &= h_i\int_0^1 P_i(\eta)\,d\eta \\ + M_i &= h_i\int_0^1 \widetilde p_i(\eta)\,d\eta \\ &= h_i\left(\frac{S_i v_0}{6} +\frac{(S_i+S_{i+1})v_m}{3} +\frac{S_{i+1}v_1}{6}\right). @@ -611,13 +636,25 @@ alias table. A common scaling of all emission values cancels from these probabilities; the implementation divides them by their maximum before forming the masses and PDFs. +Conditional on choosing interval :math:`i`, the normalized PDF with respect to +:math:`d\eta` is + +.. math:: + + p(\eta\mid i)=\frac{h_i\widetilde p_i(\eta)}{M_i}, + \qquad 0\le\eta\le1. + +The factor :math:`h_i` accounts for the change from :math:`\rho` to +:math:`\eta`; it is constant within the selected interval and cancels from +the rejection-sampling acceptance ratio. + Within the selected interval, OpenMC uses uniform proposals :math:`\eta=\xi_1` with rejection sampling. An upper bound follows by writing the cubic in the Bernstein basis: .. math:: - P_i(\eta) = \sum_{j=0}^{3}\beta_j + \widetilde p_i(\eta) = \sum_{j=0}^{3}\beta_j {3\choose j}\eta^j(1-\eta)^{3-j}, where @@ -632,13 +669,14 @@ where \end{aligned} The Bernstein basis functions are nonnegative and sum to one, so -:math:`B_i^{r}=\max_j\beta_j` bounds :math:`P_i` throughout the interval. +:math:`B_i^{\mathrm{rad}}=\max_j\beta_j` bounds :math:`\widetilde p_i` +throughout the interval. The implementation adds a small relative margin for floating-point roundoff and stores the coefficients divided by that bound. A proposal is accepted when .. math:: - \xi_2 < \frac{P_i(\eta)}{B_i^{r}}. + \xi_2 < \frac{\widetilde p_i(\eta)}{B_i^{\mathrm{rad}}}. On rejection, only :math:`\eta` is redrawn; the selected interval is retained. Selecting a new interval after a rejection would change its probability in @@ -651,30 +689,41 @@ proportion to its acceptance rate. After acceptance, Conditional Angular Sampling ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -With :math:`\rho` fixed, the target PDF is +With :math:`\rho` fixed, the target PDF over the full torus, +:math:`\theta,\zeta\in[0,2\pi)`, is .. math:: p(\theta,\zeta\mid\rho) - = \frac{f_i(\eta,\theta,\zeta)}{V'_i(\eta)}. + = \frac{J_i(\eta,\theta,\zeta)}{V'_i(\eta)}. The angular proposal uses a 32 by 32 partition of one field period, with a separate bound for each cell that varies with radius. To describe the cell -geometry, introduce the field-period angle -:math:`\varphi=N_{\mathrm{fp}}\zeta` and write +geometry, use a toroidal angle :math:`\zeta_0` restricted to the first field +period and its rescaled counterpart :math:`\varphi`: .. math:: - F_i(\eta,\theta,\varphi) - = f_i(\eta,\theta,\varphi/N_{\mathrm{fp}}). + \zeta_0\in[0,2\pi/N_{\mathrm{fp}}), \qquad + \varphi=N_{\mathrm{fp}}\zeta_0\in[0,2\pi). + +The geometric volume factor evaluated in these coordinates is + +.. math:: + + \widetilde J_i(\eta,\theta,\varphi) + = J_i(\eta,\theta,\zeta_0) + = J_i(\eta,\theta,\varphi/N_{\mathrm{fp}}). Both :math:`\theta` and :math:`\varphi` range over :math:`[0,2\pi)`, so each cell has angular half-width :math:`d=\pi/32` in both coordinates. This cell grid is separate from the finer quadrature grid used during initialization. The quadrature grid resolves the density's Fourier coefficients; the cell grid -defines regions for proposing samples. Integrating :math:`F_i` over the -rescaled square gives :math:`V'_i(\eta)`: the factor -:math:`d\varphi=N_{\mathrm{fp}}\,d\zeta` accounts for all field periods. +defines regions for proposing samples. Integrating :math:`\widetilde J_i` +over the rescaled square gives :math:`V'_i(\eta)`: the factor +:math:`d\varphi=N_{\mathrm{fp}}\,d\zeta_0` accounts for all field periods. +The normalized angular PDF in the rescaled square is therefore +:math:`p(\theta,\varphi\mid\rho)=\widetilde J_i/V'_i`. Within a fixed radial interval, express the quadratic dependence on :math:`\eta` in the Bernstein basis: @@ -683,7 +732,7 @@ Within a fixed radial interval, express the quadratic dependence on :label: stellarator-source-angular-bernstein \begin{aligned} - F_i(\eta,\theta,\varphi) + \widetilde J_i(\eta,\theta,\varphi) &= \sum_{j=0}^2 b_j(\eta)G_j(\theta,\varphi), \\ b_0(\eta)&=(1-\eta)^2, \qquad b_1(\eta)=2\eta(1-\eta), \qquad b_2(\eta)=\eta^2. @@ -695,10 +744,10 @@ values follow from the radial evaluations already available during setup: .. math:: \begin{aligned} - G_0 &= F_i(0,\theta,\varphi), \\ - G_1 &= 2F_i(1/2,\theta,\varphi) - -\tfrac12[F_i(0,\theta,\varphi)+F_i(1,\theta,\varphi)], \\ - G_2 &= F_i(1,\theta,\varphi). + G_0 &= \widetilde J_i(0,\theta,\varphi), \\ + G_1 &= 2\widetilde J_i(1/2,\theta,\varphi) + -\tfrac12[\widetilde J_i(0,\theta,\varphi)+\widetilde J_i(1,\theta,\varphi)], \\ + G_2 &= \widetilde J_i(1,\theta,\varphi). \end{aligned} This is an exact change of polynomial basis. The useful property is that @@ -710,23 +759,23 @@ then :label: stellarator-source-cell-bound B_c(\eta)=\sum_{j=0}^2 b_j(\eta)U_{jc} - \;\ge\; F_i(\eta,\theta,\varphi) + \;\ge\; \widetilde J_i(\eta,\theta,\varphi) \quad\text{within cell }c. Thus three bounds computed once per cell produce a bound at every radius in the interval. For a valid geometry, the control functions :math:`G_j` are also nonnegative. This follows from the factorization -:math:`F_i=R\widetilde\tau`, where :math:`\widetilde\tau=\sigma\tau` is the -nonnegative poloidal Jacobian. Both factors are linear in :math:`\eta`. -Writing their endpoint values as :math:`R_0,R_1` and -:math:`\widetilde\tau_0,\widetilde\tau_1` gives +:math:`\widetilde J_i=R\widetilde\tau`, where +:math:`\widetilde\tau=\sigma\tau` is the nonnegative poloidal Jacobian. +Both factors are linear in :math:`\eta`. Writing their endpoint values as +:math:`R_i,R_{i+1}` and :math:`\widetilde\tau_i,\widetilde\tau_{i+1}` gives .. math:: \begin{aligned} - G_0 &= R_0\widetilde\tau_0, \\ - G_1 &= \tfrac12(R_0\widetilde\tau_1+R_1\widetilde\tau_0), \\ - G_2 &= R_1\widetilde\tau_1. + G_0 &= R_i\widetilde\tau_i, \\ + G_1 &= \tfrac12(R_i\widetilde\tau_{i+1}+R_{i+1}\widetilde\tau_i), \\ + G_2 &= R_{i+1}\widetilde\tau_{i+1}. \end{aligned} Each expression is nonnegative because the endpoint radii are positive and @@ -901,11 +950,12 @@ when .. math:: - \xi B_c(\eta) < F_i(\eta,\theta,\varphi). + \xi B_c(\eta) < \widetilde J_i(\eta,\theta,\varphi). All cells have the same area, so the proposal density within a cell is proportional to :math:`B_c`. Multiplying it by the acceptance probability -:math:`F_i/B_c` gives a density proportional to :math:`F_i`, as required. +:math:`\widetilde J_i/B_c` gives a density proportional to +:math:`\widetilde J_i`, as required. On rejection, both the control index and the cell are redrawn while retaining the radial interval and :math:`\rho`. @@ -914,7 +964,8 @@ After acceptance, an independently and uniformly chosen field-period index .. math:: - \zeta=\frac{\varphi+2\pi\ell}{N_{\mathrm{fp}}}. + \zeta=\zeta_0+\frac{2\pi\ell}{N_{\mathrm{fp}}} + =\frac{\varphi+2\pi\ell}{N_{\mathrm{fp}}}. The sampled :math:`R` and :math:`Z` apply in every field period. The cell partition affects rejection efficiency while the accepted points follow the @@ -970,16 +1021,18 @@ distribution .. math:: - p_{\mathcal{A}}(X) = - \frac{p(X)\,\mathbf{1}_{\mathcal{A}}(X)} - {\int_{\mathcal{A}}p(Y)\,dY}, - -where :math:`X` denotes the sampled phase-space variables. The base source -wrapper repeats sampling until a candidate passes, or, with the ``kill`` -strategy, returns a failing candidate with weight zero. Independent sources -handle constraints internally: they resample positions for spatial constraints -and energies for energy constraints and the available energy range. Their time -samples are drawn directly. Constraints on an element source in a + p_{\mathcal{A}}(\boldsymbol{x}) = + \frac{p(\boldsymbol{x})\,\mathbf{1}_{\mathcal{A}}(\boldsymbol{x})} + {\int_{\mathcal{A}}p(\boldsymbol{y})\,d\boldsymbol{y}}, + +where :math:`\boldsymbol{x}=(\mathbf{r},\mathbf{\Omega},E,t)` denotes a +phase-space point and :math:`\mathbf{1}_{\mathcal A}` is one inside the +acceptance region and zero outside it. The base source wrapper repeats +sampling until a candidate passes, or, with the ``kill`` strategy, returns a +failing candidate with weight zero. Independent sources handle constraints +internally: they resample positions for spatial constraints and energies for +energy constraints and the available energy range. Their time samples are drawn +directly. Constraints on an element source in a :class:`openmc.MeshSource` are applied after selecting that element; constraints on the mesh source itself are checked by its outer wrapper. Similarly, selecting among multiple top-level sources precedes each source's constraint handling. @@ -993,12 +1046,13 @@ reasonable. Source biasing draws from a proposal distribution :math:`g` in place of the desired density :math:`p`, with the correction -:math:`w\leftarrow w\,p(X)/g(X)`. The proposal must cover every region where -the desired distribution contributes. When multiple independent factors are -biased, their corrections multiply. Mesh and point-cloud distributions can -similarly change their discrete selection probabilities and apply the ratio -of the original to the biased normalized probability. The derivation and -available biasing mechanisms are described in :ref:`methods_source_biasing`. +:math:`w\leftarrow w\,p(\boldsymbol{x})/g(\boldsymbol{x})`. The proposal must +cover every region where the desired distribution contributes. When multiple +independent factors are biased, their corrections multiply. Mesh and point-cloud +distributions can similarly change their discrete selection probabilities and +apply the ratio of the original to the biased normalized probability. The +derivation and available biasing mechanisms are described in +:ref:`methods_source_biasing`. .. [TrefethenWeideman] Lloyd N. Trefethen and J. A. C. Weideman, `The Exponentially Convergent Trapezoidal Rule From 1649ea20d5439053a4e3120676f2f1864cf80c43 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 24 Oct 2026 12:42:58 -0500 Subject: [PATCH 13/15] Clearly define joint PDF --- docs/source/methods/source_sampling.rst | 85 +++++++++++++++++-------- 1 file changed, 59 insertions(+), 26 deletions(-) diff --git a/docs/source/methods/source_sampling.rst b/docs/source/methods/source_sampling.rst index a1c0ceb0904..b0bac6855b9 100644 --- a/docs/source/methods/source_sampling.rst +++ b/docs/source/methods/source_sampling.rst @@ -449,8 +449,7 @@ cylindrical coordinates gives Here subscripts denote partial derivatives. The derivatives :math:`R_\zeta` and :math:`Z_\zeta` drop out of the determinant because the cylindrical azimuth -is exactly :math:`\zeta`. Thus the joint coordinate PDF is -:math:`p(\rho,\theta,\zeta)\propto S(\rho)R|\tau|`. +is exactly :math:`\zeta`. The geometry must have positive :math:`R` and nested flux surfaces. The sign of :math:`\tau` depends on the convention for measuring the poloidal angle: @@ -464,6 +463,46 @@ Jacobian vanishes. OpenMC determines this sign during initialization and sets :math:`J=\sigma R\tau=R|\tau|`. Initialization checks for nonpositive radii and significant changes in the sign of :math:`\tau` on a setup grid. +The role of :math:`J` is to convert a coordinate element into physical volume: + +.. math:: + + dV=J(\rho,\theta,\zeta)\,d\rho\,d\theta\,d\zeta. + +Multiplying by the emission density :math:`S(\rho)` gives the emission rate +in that element. Define the differential volume and total emission rate as + +.. math:: + + V'(\rho)=\int_0^{2\pi}\!\int_0^{2\pi} + J(\rho,\theta,\zeta)\,d\theta\,d\zeta, + \qquad Q=\int S(u)V'(u)\,du, + +where the radial integral covers the tabulated radial range. The normalized +joint PDF with respect to :math:`d\rho\,d\theta\,d\zeta` is therefore + +.. math:: + :label: stellarator-source-joint + + p(\rho,\theta,\zeta) + =\frac{S(\rho)J(\rho,\theta,\zeta)}{Q}. + +Integrating out the angles gives the radial marginal. Dividing the joint PDF +by that marginal gives the conditional angular PDF: + +.. math:: + :label: stellarator-source-factorization + + p_\rho(\rho)=\frac{S(\rho)V'(\rho)}{Q}, \qquad + p(\theta,\zeta\mid\rho) + =\frac{J(\rho,\theta,\zeta)}{V'(\rho)}. + +Thus :math:`p(\rho,\theta,\zeta)=p_\rho(\rho)p(\theta,\zeta\mid\rho)`. +OpenMC first samples :math:`\rho` from the radial marginal, then samples the +angles at that radius. The emission profile affects the radial distribution +and cancels from the conditional angular PDF. The following sections derive +the interval polynomials and angular bounds used to sample these two factors. + Interpolation Within a Radial Interval ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -500,15 +539,14 @@ factor explicitly as J_i(\eta,\theta,\zeta)=J(\rho_i+h_i\eta,\theta,\zeta). -The differential volume is the angular integral of this factor. Denote its -value at :math:`\rho=\rho_i+h_i\eta` by :math:`V'_i(\eta)`: +The differential volume defined above therefore has the interval form .. math:: :label: stellarator-source-volume-polynomial \begin{aligned} V'_i(\eta) - &\equiv \left.\frac{dV}{d\rho}\right|_{\rho=\rho_i+h_i\eta} \\ + &\equiv V'(\rho_i+h_i\eta) \\ &= \int_0^{2\pi}\!\int_0^{2\pi} J_i(\eta,\theta,\zeta)\,d\theta\,d\zeta \\ &= A_0+A_1\eta+A_2\eta^2. @@ -601,9 +639,11 @@ angular bounds needed for sampling. Selecting and Sampling a Radial Interval ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Multiplying the linear emission profile by -:eq:`stellarator-source-volume-polynomial` yields an unnormalized cubic radial -density with respect to :math:`d\rho`, expressed using the local coordinate: +To sample the radial marginal in :eq:`stellarator-source-factorization`, +multiply the linear emission profile by +:eq:`stellarator-source-volume-polynomial`. This gives an unnormalized cubic +radial density with respect to :math:`d\rho`, expressed using the local +coordinate: .. math:: :label: stellarator-source-radial-polynomial @@ -630,11 +670,11 @@ of variable :math:`d\rho=h_i\,d\eta`: +\frac{S_{i+1}v_1}{6}\right). \end{aligned} -The last expression is Simpson's rule, which is exact for a cubic. OpenMC -selects interval :math:`i` with probability :math:`M_i/\sum_j M_j` using an -alias table. A common scaling of all emission values cancels from these -probabilities; the implementation divides them by their maximum before forming -the masses and PDFs. +The last expression is Simpson's rule, which is exact for a cubic. The masses +sum to the total emission rate :math:`Q=\sum_i M_i`. OpenMC selects interval +:math:`i` with probability :math:`M_i/Q` using an alias table. A common scaling +of all emission values cancels from these probabilities; the implementation +divides them by their maximum before forming the masses and PDFs. Conditional on choosing interval :math:`i`, the normalized PDF with respect to :math:`d\eta` is @@ -681,25 +721,18 @@ and stores the coefficients divided by that bound. A proposal is accepted when On rejection, only :math:`\eta` is redrawn; the selected interval is retained. Selecting a new interval after a rejection would change its probability in proportion to its acceptance rate. After acceptance, -:math:`\rho=\rho_i+h_i\eta` has the required marginal -:math:`p_\rho\propto S V'`. +:math:`\rho=\rho_i+h_i\eta` has the required radial marginal. .. _methods_stellarator_angular_sampling: Conditional Angular Sampling ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -With :math:`\rho` fixed, the target PDF over the full torus, -:math:`\theta,\zeta\in[0,2\pi)`, is - -.. math:: - - p(\theta,\zeta\mid\rho) - = \frac{J_i(\eta,\theta,\zeta)}{V'_i(\eta)}. - -The angular proposal uses a 32 by 32 partition of one field period, with a -separate bound for each cell that varies with radius. To describe the cell -geometry, use a toroidal angle :math:`\zeta_0` restricted to the first field +With :math:`\rho` fixed, OpenMC samples the conditional angular PDF in +:eq:`stellarator-source-factorization`. The angular proposal uses a 32 by 32 +partition of one field period, with a separate bound for each cell that varies +with radius. To describe the cell geometry, use a toroidal angle +:math:`\zeta_0` restricted to the first field period and its rescaled counterpart :math:`\varphi`: .. math:: From cfc0f14517e81622b8e5d125f2a749d97b55748f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 24 Oct 2026 13:12:39 -0500 Subject: [PATCH 14/15] Fix consistency of comments --- docs/source/methods/source_sampling.rst | 7 +++---- include/openmc/source.h | 20 +++++++++++++------- openmc/source.py | 9 +++++++-- src/source.cpp | 1 - 4 files changed, 23 insertions(+), 14 deletions(-) diff --git a/docs/source/methods/source_sampling.rst b/docs/source/methods/source_sampling.rst index b0bac6855b9..68fc6c9e403 100644 --- a/docs/source/methods/source_sampling.rst +++ b/docs/source/methods/source_sampling.rst @@ -12,15 +12,14 @@ eigenvalue calculation, the external source supplies the initial population; subsequent generations use the fission bank as described in :ref:`method-successive-generations`. -This chapter describes the algorithms represented by the subclasses of +This section describes the algorithms represented by the subclasses of :class:`openmc.SourceBase`. Input examples are given in :ref:`usersguide_source`. Unless otherwise stated, :math:`\xi` denotes an independent uniform random number on :math:`[0,1)`, generated as described in :ref:`methods_random_numbers`. Normalized probability densities are denoted by :math:`p` or :math:`f`, cumulative distribution functions by :math:`F`, and -discrete selection probabilities by :math:`P`. A tilde on :math:`p` indicates -an unnormalized density. Geometric volume factors are denoted by :math:`J`. -Densities are defined with respect to the coordinates explicitly indicated. +discrete selection probabilities by :math:`P`. A tilde on :math:`p` indicates an +unnormalized density. .. _methods_source_selection: diff --git a/include/openmc/source.h b/include/openmc/source.h index c35a186b264..f84f08dd38a 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -433,15 +433,20 @@ class TokamakSource : public Source { //! by Fourier quadrature. Each radial bin is selected using its analytically //! integrated cubic PDF, then rho is sampled within the bin by rejection. //! 2. Sample (theta, zeta) from the conditional p(theta, zeta | rho) -//! ~ R |tau| by rejection against a precomputed per-radial-bin majorant. -//! 3. Sample energy and time from user-provided distribution(s). -//! 4. Sample an isotropic direction and map to Cartesian coordinates via +//! ~ R |tau| using 32 x 32 angular cells covering one field period. Select +//! a cell using a quadratic Bernstein mixture, propose angles uniformly +//! within it, and reject against its radius-dependent bound. Then choose +//! a field period uniformly. +//! 3. Map to Cartesian coordinates via //! x = R cos(zeta), y = R sin(zeta), z = Z. +//! 4. Sample an isotropic direction, then energy and time from user-provided +//! distribution(s). //! //! Fourier coefficients are interpolated linearly in rho between the surfaces -//! provided on the rho grid; within a radial bin the density R*tau is then -//! exactly quadratic in the interpolation parameter, which allows a rigorous -//! per-bin rejection majorant to be constructed from three evaluations. +//! provided on the rho grid. Within each radial bin, R*|tau| is quadratic in +//! the interpolation parameter and has three angular Bernstein control +//! functions. Taylor bounds on their Fourier series give precomputed cell +//! bounds, combined with the Bernstein weights at the sampled radius. //============================================================================== class StellaratorSource : public Source { @@ -458,7 +463,8 @@ class StellaratorSource : public Source { //========================================================================== // Private methods - //! Precompute the radial CDF and per-bin rejection majorants + //! Precompute cubic radial PDFs, bin masses, and angular cell bounds and + //! alias tables for the quadratic Bernstein mixture void precompute_sampling_distributions(); //! Evaluate the (signed) joint density f = R * tau at a point diff --git a/openmc/source.py b/openmc/source.py index d13957133e1..c36540a33f0 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -1367,8 +1367,13 @@ class StellaratorSource(SourceBase): volume for the interpolated geometry) by integrating a cubic PDF in each radial interval and using rejection sampling within the selected interval. The two angles are then sampled from the conditional - :math:`p(\theta,\zeta|\rho) \propto R|\tau|` by rejection against a - precomputed per-radial-bin majorant. + :math:`p(\theta,\zeta|\rho) \propto R|\tau|` using a 32 by 32 angular + cell partition of one field period. Within each radial interval, a + quadratic Bernstein mixture selects cells in proportion to their + radius-dependent bounds. Angles are proposed uniformly within the selected + cell and accepted by rejection sampling against its bound. The bounds are + precomputed using Taylor bounds on the Fourier series of the three + Bernstein control functions. A field period is then chosen uniformly. .. versionadded:: 0.16.1 diff --git a/src/source.cpp b/src/source.cpp index 55b79c4c438..da4bbab140d 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -14,7 +14,6 @@ #include // for dlopen, dlsym, dlclose, dlerror #endif -#include "openmc/tensor.h" #include #include "openmc/bank.h" From 17765c1e4301f5a10906b0b7eb5093b0cd89ed1e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 24 Oct 2026 13:16:04 -0500 Subject: [PATCH 15/15] Revert changes in settings.rst --- docs/source/usersguide/settings.rst | 26 -------------------------- 1 file changed, 26 deletions(-) diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 16e45cf579b..c793f815955 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -167,9 +167,6 @@ least 10,000 particles per generation to avoid any bias in the estimate of External Source Distributions ----------------------------- -The sampling algorithms and mathematical derivations for the source classes -are described in :ref:`methods_source_sampling`. - External source distributions can be specified through the :attr:`Settings.source` attribute. If you have a single external source, you can create an instance of any of the subclasses of :class:`openmc.SourceBase` @@ -396,9 +393,6 @@ A surface source file is read back like any other source file, with the Tokamak Plasma Sources ---------------------- -The volume Jacobian, radial distribution, and conditional poloidal sampling -are derived in :ref:`methods_tokamak_source`. - 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 @@ -448,9 +442,6 @@ are born at :math:`t=0`. The toroidal extent can be restricted with Stellarator Plasma Sources -------------------------- -The Fourier volume Jacobian and the radial and angular rejection algorithms -are derived in :ref:`methods_stellarator_source`. - For fully three-dimensional plasmas, the :class:`openmc.StellaratorSource` class samples neutron positions from the flux-surface Fourier representation shared by the `VMEC `_ and `DESC @@ -485,23 +476,6 @@ radial proposals are resampled within the same interval. The two angles are sampled from the conditional distribution :math:`p(\theta,\zeta|\rho) \propto R\,|\tau|` (with :math:`\tau` the poloidal-plane Jacobian) by rejection. -Angular proposals use a 32 by 32 partition of one field period in the two -spatial angular coordinates. Each cell has a bound that varies quadratically -across the radial interval. OpenMC selects -cells in proportion to these bounds, samples uniformly inside the chosen cell, -and accepts with probability equal to the density divided by the cell bound. -Rejected angular proposals retain the radius and select a new cell. A uniformly -chosen field period maps the accepted point to the full torus. The construction -of the bounds and the Bernstein mixture used to select cells are derived in -:ref:`methods_stellarator_angular_sampling`. - -The bounds are constructed during C++ initialization from the density's Fourier -coefficients, using a Taylor bound on each cell. This requires no additional -dependencies or Python preprocessing. Sampling checks for a density exceeding -its bound and reports an error if one is encountered. Bounds and alias tables -require approximately 72 KiB per radial interval on a 64-bit platform, in -addition to the geometry coefficients and temporary initialization arrays. - Sources are most conveniently created directly from equilibrium code output:: # From a VMEC wout NetCDF file