Feasibility of streamline upwind Petrov-Galerkin angular stabilization of the linear Boltzmann transport equation with magnetic fields

Radiation transport is modelled via the continuous slowing down (CSD) LBTE with magnetic field term, as was derived by St-Aubin et al (St-Aubin et al 2015, 2016). Each particle is considered based on its spatial location $\vec{r}$ and velocity $\vec{v}$. The LBTE then models the angular flux $\psi (\vec{r},\vec{{\rm{\Omega }}},E)$ parameterized by space, direction of propagation $\vec{{\rm{\Omega }}}=\vec{v}/| v|$, and energy (E). It is an integro-partial differential equation based on particle conservation in six-dimensional phase space. The magnetic field impacts charged particles due to the Lorentz force, which is modelled via an angular streaming term (St-Aubin et al 2015, 2016).

In external beam mega-voltage (MV) photon radiotherapy, electrons are generated from photon-electron interactions, and subsequently by electron-electron interactions. Thus, in order to fully characterize the system, a coupled system of equations is required. In setting up the coupled system, two assumptions are made: first, bremsstrahlung interactions are neglected, a reasonable assumption in the MV therapeutic energy range (Bush et al 2011), and second, positrons that are produced do not annihilate, and are thus treated as electrons (St-Aubin et al 2015). These assumptions allow the photon equation to be solved, followed by the electron equation. The resulting coupled system is given by

$$\begin{eqnarray}&&\left[\vec{{\rm{\Omega }}}\cdot {\rm{\nabla }}+{\sigma }_{\gamma }(\vec{r},E)\right]{\psi }_{\gamma }(\vec{r},\vec{{\rm{\Omega }}},E)={Q}^{\gamma \gamma }+{S}_{\gamma },\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{l}\left[\vec{{\rm{\Omega }}}\cdot {\rm{\nabla }}+\displaystyle \frac{q}{| \vec{p}| }\vec{\tau }\cdot {{\rm{\nabla }}}_{{\rm{\Omega }}}+{\sigma }_{e}(\vec{r},E)\right]{\psi }_{e}(\vec{r},\vec{{\rm{\Omega }}},E)\\ -\,\displaystyle \frac{\partial }{\partial E}({\beta }_{r}(\vec{r},E){\psi }_{e}(\vec{r},E,\vec{{\rm{\Omega }}}))={Q}^{{ee}}+{Q}^{\gamma e}+{S}_{e},\end{array}\end{eqnarray} \tag{ 2 }$$

where γ indicates photon flux, and e indicates electron flux. The charge of the particles being modelled is given by q, $\vec{p}$ is the momentum vector, σ_x is the total interaction cross-section for particle x, β_r is the restricted mass stopping power, and

$$\begin{eqnarray}&&\begin{array}{l}\vec{\tau }(\vec{B},\vec{{\rm{\Omega }}})=\displaystyle \frac{1}{\sqrt{1-{\cos }^{2}\theta }}\\ \times \,\left[{\left(\vec{{\rm{\Omega }}}\times (\vec{{\rm{\Omega }}}\times \vec{B})\right)}_{z}\hat{\phi }-{\left(\vec{{\rm{\Omega }}}\times \vec{B}\right)}_{z}\hat{\theta }\right],\end{array}\end{eqnarray} \tag{ 3 }$$

where θ and ϕ are the azimuthal and polar angles, respectively, and the magnetic field is given by $\vec{B}$. The right hand side terms of equations (1) and (2) represent the external source term S_x and the scattering source Q^xy , representing particles of type y generated due to scatter with particles of type x (St-Aubin et al 2016), given by

$$\begin{eqnarray}&&\begin{array}{l}{Q}^{{xy}}={\displaystyle \int }_{{\rm{\Omega }}}d\vec{{\rm{\Omega }}}^{\prime} {\displaystyle \int }_{0}^{\infty }{dE}^{\prime} {\sigma }_{s}^{{xy}}(\vec{r},E^{\prime} \to E,\vec{{\rm{\Omega }}}\cdot \vec{{\rm{\Omega }}}^{\prime} )\\ \times \,{\psi }_{x}(\vec{r},E^{\prime} ,\vec{{\rm{\Omega }}}^{\prime} ),\end{array}\end{eqnarray} \tag{ 4 }$$

where ${\sigma }_{s}^{{xy}}$ gives the macroscopic differential cross section (Gifford et al 2006). In this work, the focus is on angular stabilization, hence the focus will now be on the electron equation, and the subscript e is omitted.

Finally, a vacuum boundary condition is used, where

$$\begin{eqnarray}&&\psi =0,\hspace{4mm}\vec{{\rm{\Omega }}}\cdot \vec{n}\lt 0,\end{eqnarray} \tag{ 5 }$$

and $\vec{n}$ is the outward-facing unit normal to the boundary. Dose $D(\vec{r})$ is calculated from the solved electron fluence by

$$\begin{eqnarray}&&D(\vec{r})={\int }_{0}^{\infty }{dE}{\int }_{{\rm{\Omega }}}d\vec{{\rm{\Omega }}}\displaystyle \frac{{\sigma }_{{ED}}(\vec{r},E)}{\rho (\vec{r})}{\psi }_{e}(\vec{r},\vec{{\rm{\Omega }}},E),\end{eqnarray} \tag{ 6 }$$

where $\rho (\vec{r})$ is the material density and ${\sigma }_{{ED}}(\vec{r},E)$ is the energy deposition cross-section (St-Aubin et al 2016).

Due to the complexity of equation (2), the solution generally must be determined numerically. A common solution method for the LBTE is a Gauss-Seidel technique known as source iteration (Miller and Lewis 1996, Larsen and Morel 2010). In the application of this technique, the LBTE with magnetic fields is reformulated after a spherical harmonic expansion as (St-Aubin et al 2015)

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal L }{\psi }^{(t+1)}(\vec{r},\vec{{\rm{\Omega }}},E)\\ =\,\displaystyle \sum _{{\ell }=0}^{{L}_{{\rm{\max }}}}\displaystyle \sum _{m=-{\ell }}^{{\ell }}{\displaystyle \int }_{0}^{\infty }{dE}^{\prime} {\sigma }_{{\ell }}(\vec{r},E^{\prime} \to E)\\ \ \ \times \,{\phi }_{{\ell }m}^{(t)}(\vec{r},E^{\prime} ){Y}_{{\ell }m}(\vec{{\rm{\Omega }}})+S(\vec{r},\vec{{\rm{\Omega }}},E),\end{array}\end{eqnarray} \tag{ 7 }$$

where the σ_ℓ are the coefficients of a Legendre expansion of the differential scattering cross section (Zelyak 2017), and ${\phi }_{{\ell }m}^{(t)}$ is given by

$$\begin{eqnarray}&&{\phi }_{{\ell }m}^{(t)}(\vec{r},E)=\int d\vec{{\rm{\Omega }}}{\psi }^{(t)}(\vec{r},E,\vec{{\rm{\Omega }}}){Y}_{{\ell }m}^{* }(\vec{{\rm{\Omega }}}).\end{eqnarray} \tag{ 8 }$$

${ \mathcal L }$ is the operator acting on ψ on the left hand side of equation (2). The unknown ψ is then solved iteratively with iteration index t until the relative error achieves the desired precision of 10⁻⁴ for this work.

The convergence properties of the source iteration method are important in understanding the computational speed of solution. Convergence properties of a stationary iterative scheme can be quantified through the spectral radius ρ, or the absolute value of the eigenvalue with largest magnitude of the iteration matrix (Larsen and Morel 2010, Isaacson and Keller 1966). When the spectral radius remains less than 1 regardless of parameter values, the scheme is considered unconditionally stable. In this work, we consider two methods for determining the spectral radius. A Fourier analysis (Larsen and Morel 2010) is first presented, where the dependence of the spectral radius on the relevant parameters is investigated. The results of the Fourier analysis are then confirmed numerically by calculating

$$\begin{eqnarray}&&\rho ={{\rm{lim}}}_{t\to \infty }\displaystyle \frac{| | {\psi }^{(t)}-{\psi }^{(t-1)}| | }{| | {\psi }^{(t-1)}-{\psi }^{(t-2)}| | },\end{eqnarray} \tag{ 9 }$$

using the simulation code developed in this work. For practical purposes, t is taken to be large enough for the value of ρ to stabilize within a 10⁻⁶ tolerance.

Zelyak et al (Zelyak et al 2018) performed a Fourier analysis for the DGFEM in order to quantify the dependence of the spectral radius on parameters such as magnetic field strength, total cross-section, and degree of anisotropy. Their analysis verified that in the case of zero magnetic field the spectral radius expression was reduced to the ratio of the isotropic scattering cross-section ${\sigma }_{0}(\vec{r})$ to the total cross-section $\sigma (\vec{r})$:

$$\begin{eqnarray}&&\rho =\displaystyle \frac{{\sigma }_{0}(\vec{r})}{\sigma (\vec{r})}.\end{eqnarray} \tag{ 10 }$$

For the discrete DGFEM LBTE with magnetic field, it was found that the spectral radius increased with magnetic field strength yet always remained less than one. However, it was also shown that as the material density decreased, the spectral radius rapidly approached unity with increasing magnetic field strength. This result shows that although the DGFEM is unconditionally stable, it suffers from some degradation of the convergence rate as the magnetic field strength increases or as the material density decreases. In addition, the angular elements must be solved in a particular order for the DGFEM upwind stabilization scheme lest the method suffer further spectral radius degradation, which poses challenges for parallelization (Zelyak et al 2018).

Each variable must be discretized in order to solve the LBTE with magnetic fields numerically. The energy variable is first discretized using the multigroup method (Miller and Lewis 1996). In the original formulation by St-Aubin et al (St-Aubin et al 2015, 2016) and Yang et al (Yang et al 2018), both space and angle were discretized using a DGFEM for stabilization. In the work of Yang et al (Yang et al 2018), the spatial elements are given by cubic voxels with eight nodes each, and the angular domain is divided into triangles conformal to the unit sphere. The unknown angular fluence is then expanded over the basis functions for each pair of spatial and angular elements:

$$\begin{eqnarray}&&\psi \approx \tilde{\psi }=\displaystyle \sum _{i=1}^{I}\displaystyle \sum _{p=1}^{P}{\psi }_{{ip}}{\lambda }_{i}(\vec{r}){\gamma }_{p}(\vec{{\rm{\Omega }}}),\end{eqnarray} \tag{ 11 }$$

where ${\lambda }_{i}(\vec{r})$ and ${\gamma }_{p}(\vec{{\rm{\Omega }}})$ are the spatial and angular basis functions, respectively. A Galerkin method is then applied discontinuously, where the system is solved over each pair of spatial and angular nodes individually. Since the solution can be discontinuous over element boundaries, the value of ψ_ip at a given node is permitted to be different for each element. Allowing for discontinuities in this manner better allows for sharp changes in the fluence to be captured, particularly at material interfaces. Excellent accuracy was achieved via the DGFEM when compared to GEANT4 Monte Carlo simulations for a variety of slab and patient geometries (Yang et al 2018). As such, the results of the alternative SUPG angular stabilization scheme considered here will be compared to the existing validated DGFEM scheme.

An SUPG method is developed for the angular domain, with the main goal of reducing the spectral radius of the source iteration scheme, and potentially providing a more obvious pathway to parallelization. This method is applied by incorporating a modification into the weighting function effectively introducing a small amount of artificial diffusion to provide stabilization in the angular variable in the presence of advection. All previous applications of the SUPG method to the LBTE have been done using the spatial advection term (Pain et al 2006, Brooks and Hughes 1982, Merton 2012, Hughes et al 1986, Kanschat 1998), and application to angular advection would be a novel advancement.

To simplify the derivation, the magnetic field orientation will be assumed to be in the z-direction. This can be done without loss of generality as any problem can be reduced to this orientation through a coordinate transformation (St-Aubin et al 2016). Then the multigroup LBTE with magnetic field term is given by

$$\begin{eqnarray}&&\begin{array}{l}\vec{{\rm{\Omega }}}\cdot {\rm{\nabla }}\psi (\vec{r},\vec{{\rm{\Omega }}})+\kappa {B}_{z}\displaystyle \frac{\partial }{\partial \varphi }\psi (\vec{r},\vec{{\rm{\Omega }}})\\ +\,{\sigma }_{T}(\vec{r})\psi (\vec{r},\vec{{\rm{\Omega }}})=Q(\vec{r},\vec{{\rm{\Omega }}})+S(\vec{r},\vec{{\rm{\Omega }}}).\end{array}\end{eqnarray} \tag{ 12 }$$

where the CSD operator is combined with the source term Q using a diamond difference discretization (Morel 1985). The weak formulation incorporating the standard SUPG modification is

$$\begin{eqnarray}&&\left\langle { \mathcal L }\psi ,{\rm{\Upsilon }}+\delta \kappa {B}_{z}\displaystyle \frac{\partial {\rm{\Upsilon }}}{\partial \varphi }\right\rangle =\left\langle Q+S,{\rm{\Upsilon }}+\delta \kappa {B}_{z}\displaystyle \frac{\partial {\rm{\Upsilon }}}{\partial \varphi }\right\rangle ,\end{eqnarray} \tag{ 13 }$$

where ${ \mathcal L }$ is the operator acting on ψ on the left hand side of equation (12), ${\rm{\Upsilon }}={\lambda }_{j}(\vec{r}){\gamma }_{q}(\vec{{\rm{\Omega }}})$ is the original unmodified weighting function used in the DGFEM, κ is the multigroup magnetic field parameter, given by (Yang et al 2018, Zelyak et al 2018, St-Aubin et al 2016)

$$\begin{eqnarray}&&\begin{array}{l}\kappa =\displaystyle \frac{{qc}}{{E}_{g}-{E}_{g+1}}\\ \times \,\mathrm{ln}\left(\displaystyle \frac{{E}_{g}+{m}_{0}{c}^{2}+\sqrt{{\left({E}_{g}+{m}_{0}{c}^{2}\right)}^{2}-{\left({m}_{0}{c}^{2}\right)}^{2}}}{{E}_{g+1}+{m}_{0}{c}^{2}+\sqrt{{\left({E}_{g}+{m}_{0}{c}^{2}\right)}^{2}-{\left({m}_{0}{c}^{2}\right)}^{2}}}\right)\end{array}\end{eqnarray} \tag{ 14 }$$

and δ is the SUPG free parameter (Kanschat 1998, Fries and Matthies 2004). The value of this free parameter must be determined.

In order to investigate the accuracy of this technique, the fluence and subsequent dose were calculated using the derived SUPG method for a 2 × 2 cm² 6 MV photon spectrum beam (St-Aubin et al 2015) originating from a point source with an SSD of 100cm in a 9.6 × 9.6 × 9.6 cm³ phantom. This field size was selected because a small field with its associated charged particle disequilibrium makes for a challenging validation case. 2 mm Cartesian hexahedral spatial voxels were used in the beam and penumbra, 4 mm voxels were used in the volume surrounding the beam, and 8 mm voxels were used at the edges of the phantom where the fluence values were almost negligible based on the accuracy results from Yang et al (Yang et al 2018). The angular domain was divided into 32 triangular angular elements conformal to the unit sphere. Linear basis functions were used in space, while quadratic basis functions were used in angle. Energy discretization was performed using 32 photon energy groups and 20 electron energy groups and the photon source spectrum presented in the work of St-Aubin et al was used (St-Aubin et al 2015). A logarithmic Multigroup energy discretization was performed such that the energy group width decreased as the electron energy decreased. This allowed for an increased energy resolution at lower energies to account for the increasing number of interactions. At the highest energy the multigroup energy width was 1.84MeV while for the lowest energy group the width was 0.04 MeV. According to equation (14), the multigroup magnetic field parameter varies with multigroup energy. Thus, when multiplied by the magnetic field the magnetic field parameter can vary from 0.2 cm⁻¹ at the highest energy to 40.5 cm⁻¹ at the lowest energy for the field strength range considered here. The material cross-sections were generated via the CEPXS software (Lorence et al 1989) with the Legendre expansion truncated at L = 5. An investigation into the impact of changing angular mesh sizes was also performed to validate the consistency of this formalism over a range of angular mesh sizes including 8, 32 and 128 angular elements conformal to the unit sphere with quadratic basis functions.

The results were compared to the previously validated DGFEM results (Yang et al 2018), as this method showed a high degree of agreement with Monte Carlo, and it is truly this method with which the SUPG scheme is being compared in this instance. Homogeneous phantoms of lung (2.6 × 10⁻¹ gcm⁻³), water (1.00 g cm⁻³) and bone (1.92 g cm⁻³) were first used along with magnetic field configurations of a 0.5 T parallel field and a 1.5 T perpendicular field. Two heterogeneous slab phantoms were then considered. The heterogeneous water-bone-lung-water phantom (WBLW) contained 4 cm of water, followed by 1.2 cm of bone, then 2.8 cm of lung, and finally 1.6 cm of water. The heterogeneous water-bone-air-water phantom (WBAW) contained 4 cm of water, followed by 1.2 cm of bone, then 2.8 cm of air (1.2 × 10⁻³ g cm⁻³), and finally 1.6 cm of water.

In order to evaluate the accuracy, both central axis depth dose curves and lateral depth dose profiles were generated for both the SUPG method and the DGFEM in angle using 32 quadratic angular elements. Three dimensional gamma analyses at both the 2%/2 mm and 3%/3 mm levels were performed using the global Van Dyk normalization, with only points having at least 10% of the maximum dose qualifying for the analysis.

The spectral radius of the SUPG method was analyzed both using the Fourier analysis method and the asymptotic behaviour of the simulation code. The adjoint SUPG operator is given by ${A}^{\dagger }=1-\delta \kappa B\tfrac{\partial }{\partial \varphi }$, and the continuous SUPG equation is given by applying this adjoint operator to equation (12):

$$\begin{eqnarray}&&{A}^{\dagger }{ \mathcal L }\psi ={A}^{\dagger }g,\end{eqnarray} \tag{ 15 }$$

where ${ \mathcal L }$ is as defined in equation (13), and $g=Q(\vec{r},\vec{{\rm{\Omega }}})+S(\vec{r},\vec{{\rm{\Omega }}})$ as given in equation (12). We first defined error functions and their Fourier expansions as

$$\begin{eqnarray}&&\begin{array}{l}\delta {\psi }^{(t+1)}(\vec{r},\vec{{\rm{\Omega }}})={\psi }^{(t+1)}(\vec{r},\vec{{\rm{\Omega }}})-{\psi }^{(t)}(\vec{r},\vec{{\rm{\Omega }}})\\ =\,\displaystyle \sum _{p}{\displaystyle \int }_{-\infty }^{\infty }d\vec{\lambda }{a}_{p}^{(t+1)}(\vec{\lambda },\theta ){e}^{i\vec{\lambda }\vec{r}}{e}^{{ip}\varphi },\end{array}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\begin{array}{l}\delta {\phi }_{{\ell }m}^{(t+1)}(\vec{r})={\phi }_{{\ell }m}^{(t+1)}(\vec{r})-{\phi }_{{\ell }m}^{(t)}(\vec{r})\\ \quad =\,{\displaystyle \int }_{-\infty }^{\infty }d\vec{\lambda }{b}_{{\ell }m}^{(t+1)}(\vec{\lambda }){e}^{i\vec{\lambda }\vec{r}}.\end{array}\end{eqnarray} \tag{ 17 }$$

Then the iterative system of equations incorporating the error functions was given by

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{p=-\infty }^{\infty }(1-\delta \kappa {Bp})[i\vec{{\rm{\Omega }}}\cdot \vec{\lambda }+i\kappa {Bp}\\ \quad +\,\sigma ]{a}_{p}^{(t+1)}(\vec{\lambda },\theta ){e}^{{ip}\varphi }\\ \ =\,\left(1-\delta \kappa B\displaystyle \frac{\partial }{\partial \varphi }\right)\displaystyle \sum _{{\ell }m}{\sigma }_{{\ell }}{Y}_{{\ell }m}(\vec{{\rm{\Omega }}}){b}_{{\ell }m}^{(t)}(\vec{\lambda }),\end{array}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{b}_{{\ell }m}^{(t)}(\vec{\lambda })={\int }_{{\rm{\Omega }}}d\vec{{\rm{\Omega }}}{Y}_{{\ell }m}^{* }(\vec{{\rm{\Omega }}})\displaystyle \sum _{m^{\prime} =\infty }^{\infty }{a}_{m^{\prime} }^{(t)}(\vec{\lambda },\theta ){e}^{{im}^{\prime} \varphi }.\end{eqnarray} \tag{ 19 }$$

equations (16) and (17) were substituted into equations (18) and (19) and simplified (Zelyak et al 2018, Zelyak 2017) to obtain the matrix equation

$$\begin{eqnarray}&&{b}_{{\ell }m}^{(t)}(\vec{\lambda })={T}_{{\ell }{\ell }^{\prime} }{b}_{{\ell }m}^{(t-1)}(\vec{\lambda }),\end{eqnarray} \tag{ 20 }$$

with matrix components given by

$$\begin{eqnarray}&&{T}_{{\ell }{\ell }^{\prime} }={\sigma }_{{\ell }}^{\prime} {C}_{{\ell }{\ell }^{\prime} }^{m}{\int }_{-1}^{1}\displaystyle \frac{d\mu }{2}\displaystyle \frac{{P}_{{\ell }m}(\mu ){P}_{{\ell }^{\prime} m}(\mu )}{i\lambda \mu +i\kappa {Bm}+\sigma },\end{eqnarray} \tag{ 21 }$$

and

$$\begin{eqnarray}&&{C}_{{\ell }{\ell }^{\prime} }^{m}=\sqrt{\displaystyle \frac{(2{\ell }+1)({\ell }-m)!}{({\ell }+m)!}}\sqrt{\displaystyle \frac{(2{\ell }^{\prime} +1)({\ell }^{\prime} -m)!}{({\ell }^{\prime} +m)!}}.\end{eqnarray} \tag{ 22 }$$

At this point, the key result was that the SUPG operator completely canceled out and thus had no effect on the spectral radius compared to the original standard Galerkin formulation. In particular, the final result of equation (21) was identical to that of Zelyak et al (Zelyak et al 2018, Zelyak 2017). Through further analysis, Zelyak et al showed that the spectral radius (absolute value of the largest eigenvalue) for the matrix system of equation (20) reduced to equation (10), and thus the magnetic field has no impact on the spectral radius for the SUPG method.

These results were confirmed numerically using equation (9). The results of the continuous derivation indicate that the spectral radius should not depend on either B or δ, and that for all values of these parameters we should see ρ = c = σ₀/σ. To investigate this, a variety of parameter values were tested, including c=0.05, 0.1, 0.2, 0.4, 0.99, δ = 0.05, 0.5, 5, and κ B = 1, 5, 10, 20, 40. It was found that in all cases, ρ = c as expected.

For the DGFEM, the spectral radius had a dependence on both the magnetic field strength and the total cross-section value which is proportional to the material density. This can be seen in figure 1, which also shows the constant value of the spectral radius for the SUPG method, calculated for δ = 0.1. Realistic values of κ B range from 0.2-40.5 for the magnetic field strengths considered here.

Zoom In Zoom Out Reset image size

In order to apply the SUPG method, the value of δ needed to be determined. An empirical investigation found that in general, both higher energies and sparser materials required a higher value of δ. This can be explained by the fact that a smaller cross-section value led to a system that was more advection-dominated, and thus required additional stabilization. Previous investigations (Fries and Matthies 2004, Brooks and Hughes 1982) determined a formula for δ which is based on the coth formula, which qualitatively fit with our empirical data. The published formula however was derived based on advection-diffusion equations, where in this application there is no explicit diffusion coefficient. As such, we defined a 'pseudo-Peclet number' used to define δ:

$$\begin{eqnarray}&&\tilde{{Pe}}=\tilde{C}\displaystyle \frac{c{\rm{\Delta }}{\rm{\Omega }}}{2{\sigma }_{T}},\end{eqnarray} \tag{ 23 }$$

where c is the angular advective velocity, given by κ B_z for the z-magnetic field orientation, and ΔΩ quantifies the angular step size via the length of the triangular bisector. In order to make $\tilde{{Pe}}$ non-dimensional, a constant $\tilde{C}$ is introduced having units of $\tfrac{1}{{\rm{rad}}}$. Then the SUPG free parameter was given by

$$\begin{eqnarray}&&\delta =\alpha \displaystyle \frac{{\rm{\Delta }}{\rm{\Omega }}}{2c}\left({\rm{\coth }}\left(\displaystyle \frac{1}{\tilde{{Pe}}}\right)-\tilde{{Pe}}\right),\end{eqnarray} \tag{ 24 }$$

where the constant of proportionality α and constant $\tilde{C}$ were used to fit the function against our empirical data. For a proof of concept investigation, αΔΩ was taken to be equal in value to $\tilde{C}{\rm{\Delta }}{\rm{\Omega }}$, which henceforth will simply be referred to collectively as C_δ . This reduces the number of parameters to fit and simplifies the optimization, though potential improvements could later be made by separating the variables. A value of C_δ was determined for each material and each magnetic field strength for a fixed value of 32 angular elements with quadratic basis functions. The effect of changing the angular element size will be discussed in section 3.4. A fit was then determined for each magnetic field configuration of C_δ versus material density, shown in figure 2. A higher value of C_δ was required for the higher magnetic field strength due to the increased angular advection. An investigation of more magnetic field strength values could possibly yield an analogous dependence on field strength but only 0.5 T and 1.5 T field strengths were investigated here.

Zoom In Zoom Out Reset image size

Simulations were performed for both the homogeneous phantoms and the slab geometry phantoms and the percentage of points passing a 2%/2 mm three dimensional gamma analysis as compared to the DGFEM results are presented in table 1, as well as those passing a 3%/3 mm gamma analysis. Only points having a dose of more than 10% of the maximum dose were included in the analysis.

The results for a 0.5 T parallel magnetic field in the WBLW phantom are given in figure 3 and table 1. The maximum percent difference for depths greater than ${d}_{{\rm{\max }}}$ was 0.27%. The same phantom was then used for a 1.5 T perpendicular magnetic field with results given in figure 4 and table 1. The maximum percent difference at depths greater than ${d}_{{\rm{\max }}}$ in this case was 14.7%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results for a 0.5 T parallel magnetic field in the WBAW phantom are given in figure 5 and table 1. The maximum percent difference at depths greater that ${d}_{{\rm{\max }}}$ for the 0.5 T parallel case was 32.9%. For the 1.5 T perpendicular field case, a maximum percent difference at depths greater than ${d}_{{\rm{\max }}}$ of 52.0% was observed. These results can be seen in figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Convergence results are given in table 2, showing the reduced number of iterations required to fully converge on the final solution for all energy groups. For the denser materials, the benefit of the SUPG method on convergence is small while for the lower density materials (e.g. air) the SUPG method results in a 40% reduction in the number of iterations for the perpendicular magnetic field case.

For the investigation above, a fixed number of 32 angular elements was used to facilitate fair comparison with the reference data, however this angular mesh can be further refined to 128 elements, or coarsened to 8 elements, representing a range of angular mesh sizes that are clinically relevant (Yang et al 2018). Based on equations (23) and (24), the values of C_δ change due the differing bisector lengths for the 8, 32 and 128 angular meshes. By applying these changes to C_δ , and the resulting changes in the free parameter δ, PDDs were simulated for 8, 32, and 128 quadratic angular meshes in the WBLW slab geometry phantom for a 0.5 T parallel magnetic field. This analysis was chosen since the presented SUPG method is most promising for this configuration. The results are shown in figure 7, alongside the zero magnetic field case for comparison. Figure 7 (b) shows a zoomed in view to better see the difference between the PDDs. As was reported previously (Yang et al 2018), there is a significant improvement in accuracy when moving from 8 angular elements to 32, however the accuracy improvement is negligible in moving from 32 to 128 elements. These results show that the presented SUPG formalism is generally applicable over the angular mesh sizes investigated here which would be the most useful clinically.

Zoom In Zoom Out Reset image size