Native parametric tokamak source

[eepeterson]

Description

This PR implements a new type of source based on a common parametric description of a tokamak plasma. A schematic of the tokamak shaping parameters can be seen in the image below where $R_0$ is the major radius, $a$ is the minor radius, $\delta$ is the triangularity, $\kappa$ is the elongation, $\Delta$ is the Shafranov shift, $\tilde{r}$ is the normalized minor radial coordinate $r/a$, and $\alpha$ is the poloidal angle measured counterclockwise from the outboard midplane.

The parameterization defines the flux surfaces with the following equations:

$$\large \begin{aligned} R\left(r,\alpha\right) &= R_0 + r\cos\left(\alpha + \delta\sin\alpha\right) +\Delta\left(1-\left(\dfrac{r}{a}\right)^2\right)\\ Z\left(r,\alpha\right) &= \kappa r\sin{\alpha}, \end{aligned}$$

The total source emission rate, $S_0$ in neutrons / second is given by the integral over the plasma volume of the source emission density $S(R,Z)$ as shown in the equation below.

$$\large S_0 = \int_VS\left(R,Z\right)RdRdZd\phi$$

To sample the plasma conformally and continuously without rejection we want to transform coordinates from $R, Z$ to $r, \alpha$ because the plasma is defined to be on $r \in [0,a] \times \alpha \in [0, 2\pi]$ and the source emission density is a flux function and therefore only a function of $r$. This means the spatial probability distribution we need to sample from is (assuming axisymmetry so that $\phi$ can be sampled uniformly):

$$\large p(\tilde{r},\alpha) \propto S(\tilde{r})R(\tilde{r}, \alpha)\mathcal{J}(\tilde{r}, \alpha)$$

where $R(\tilde{r}, \alpha)$ and $\mathcal{J}(\tilde{r}, \alpha)$ are the major radius and jacobian of the transformation respectively. I won't work through the remainder of the math here, but the algorithm to sample the neutron source site proceeds as follows

1. Sample $\tilde{r}$ from the marginal distribution $p(\tilde{r})$
2. Sample $\alpha$ from the conditional distribution $p(\alpha | \tilde{r})$
3. Sample $\phi$ uniformly over the toroidal extent of the plasma
4. Sample the neutron energy
5. Sample the angle isotropically
6. Convert $(\tilde{r}, \alpha, \phi)$ coordinates to $(x, y, z)$ and return the source site.

Fixes # (issue)

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[paulromano]

Great work on this new source class @eepeterson! I know there are a lot of people in the fusion world who will be excited to use this. See my comment below about the radial distribution sampling; once we resolve that, this should be good to merge.

[paulromano]

@eepeterson I've just made my updated. I've implemented a grid refinement on the radial sampling as proposed. I've also refactored how the radial and poloidal distributions are stored so that we can take advantage of the Tabular class. One other thing that needed fixing (outside the scope of your changes) to make it work well was that Tabular used a linear search over CDF points, which I've changed to use a binary search.

It may be worth re-running some of the test cases from your manuscript just to make sure my changes haven't messed anything up, but I'll leave it up to you! As far as I'm concerned, this PR is good to go.