Angular Quadrature Research Articles

Hybrid methods for radiation transport using diagonally implicit Runge–Kutta and space–time discontinuous Galerkin time integration

In this work, we describe extensions of a hybrid method for time-dependent linear, kinetic radiation transport problems to high-order time integration schemes of the diagonally-implicit Runge–Kutta (DIRK) and space–time discontinuous Galerkin (STDG) types. The hybrid methods are constructed by splitting the radiation flux into collided and uncollided components to which low- and high-resolution discrete ordinates approximations are applied. The efficiency of hybrid methods constructed using DIRK, STDG, integral deferred correction (IDC), and implicit Euler schemes is compared using a test problem in one-dimensional slab geometry containing material discontinuities. It is observed that (i) higher-order methods are more efficient than the implicit Euler method, often by an order of magnitude or more; (ii) third-order methods yield solutions of a given error in roughly half the time of the second-order method of the same type; and (iii) for a given order of accuracy it is found that the most efficient class of time integration scheme is STDG, followed by DIRK, with IDC the least efficient, for the test problem considered. Two test problems in two-dimensional xy-geometry are used to compare the computational efficiency of hybrid and standard discrete ordinates methods constructed with DIRK and STDG integrators. We observe that replacing a standard discrete ordinates method using an angular quadrature of order N with a hybrid discrete ordinates method using angular quadratures of order 2N and N/2 for the uncollided and collided fluxes, respectively, usually reduces overall solution time by a factor of 2 or more while simultaneously reducing the resulting solution error by a factor of 2 or more for the test problems considered.

An efficient pseudo-spectral method for the description of atomic electronic wave functions – Application to the hydrogen atom in a uniform magnetic field

The mapping of an electronic state on a real-space support lattice may offer advantages over a basis set ansatz in cases where there are linear dependences due to basis set overcompleteness or when strong internal or external fields are present. Such discretization methods are also of interest because they allow for the convenient numerical integration of matrix elements of local operators. We have developed a pseudo-spectral approach to the numerical solution of the time-dependent and time-independent Schrödinger equations describing electronic motion in atoms and atomic ions in terms of a spherical coordinate system. A key feature of this scheme is the construction of a Variational Basis Representation (VBR) for the non-local component and of a Generalized Finite Basis Representation (GFBR) for the local component of the electronic Hamiltonian operator. Radial Hamiltonian eigenfunctions χnl;βr of the H atom-like system and spherical harmonics form the basis set. Two special cases of this approach are explored: In one case, the functions of the field-free H atom are used as the elements of the basis set, and in the second case, each radial basis function has been obtained by variationally optimizing a shielding parameter β to yield a minimum energy for a particular eigenstate of the H atom in a uniform magnetic field.We derive a new quadrature rule of nearly Gaussian accuracy for the computation of matrix elements of local operators between radial basis functions and perform numerical evaluation of local operator matrix elements involving spherical harmonics. With this combination of radial and angular quadrature prescriptions we ensure to a good approximation the discrete orthogonality of Hamiltonian eigenfunctions of H atom-like systems for summation over the grid points. We further show that sets of χnl;βr functions are linearly independent, whereas sets of the polar-angle-dependent components of the spherical harmonics, i.e., the associated Legendre functions, are not and provide a physical interpretation of this mathematical observation.The pseudo-spectral approach presented here is applied to two model systems: the field-free H atom and the H atom in a uniform magnetic field. The results demonstrate the potential of this method for the description of challenging systems such as highly charged atomic ions.

The adaptive collision source method for discrete ordinates radiation transport

A novel collision source method has been developed to solve the Linear Boltzmann Equation (LBE) more efficiently by adaptation of the angular quadrature order. The angular adaptation method is unique in that the flux from each scattering source iteration is obtained, with potentially a different quadrature order used for each. Traditionally, the flux from every iteration is combined, with the same quadrature applied to the combined flux. Since the scattering process tends to distribute the radiation more evenly over angles (i.e., make it more isotropic), the quadrature requirements generally decrease with each iteration. This method allows for an optimal use of processing power, by using a high order quadrature for the first iterations that need it, before shifting to lower order quadratures for the remaining iterations. This is essentially an extension of the first collision source method, and is referred to as the adaptive collision source (ACS) method. The ACS methodology has been implemented in the 3-D, parallel, multigroup discrete ordinates code TITAN. This code was tested on a several simple and complex fixed-source problems. The ACS implementation in TITAN has shown a reduction in computation time by a factor of 1.5–4 on the fixed-source test problems, for the same desired level of accuracy, as compared to the standard TITAN code.

Discrete Ordinates Quadratures Based on Linear and Quadratic Discontinuous Finite Elements over Spherical Quadrilaterals

We present a new family of discrete ordinates (Sn) angular quadratures based on discontinuous finite elements (DFEMs) in angle. The angular domain is divided into spherical quadrilaterals (SQs) on the unit sphere surface. Linear discontinuous finite element (LDFE) and quadratic discontinuous finite element (QDFE) basis functions in the direction cosines are defined over each SQ, producing LDFE-SQ and QDFE-SQ angular quadratures, respectively. The new angular quadratures demonstrate more uniform direction and weight distributions than previous DFEM-based angular quadratures, local refinement capability, strictly positive weights, generation to large numbers of directions, and fourth-order accurate high-degree spherical harmonics (SH) integration. Results suggest that particle-conservation errors due to inexact high-degree SH integration rapidly diminish with quadrature refinement and tend to be orders of magnitude smaller than other discretization errors affecting the solution. Results also demonstrate that the performance of the new angular quadratures without local refinement is on par with or better than that of traditional angular quadratures for various radiation transport problems. The performance of the new angular quadratures can be further improved by using local refinement, especially within an adaptive Sn algorithm.

Evaluation of a deterministic grid-based Boltzmann solver (GBBS) for voxel-level absorbed dose calculations in nuclear medicine

To evaluate the 3D Grid-based Boltzmann Solver (GBBS) code ATTILA® for coupled electron and photon transport in the nuclear medicine energy regime for electron (beta, Auger and internal conversion electrons) and photon (gamma, x-ray) sources. Codes rewritten based on ATTILA are used clinically for both high-energy photon teletherapy and 192Ir sealed source brachytherapy; little information exists for using the GBBS to calculate voxel-level absorbed doses in nuclear medicine.We compared DOSXYZnrc Monte Carlo (MC) with published voxel-S-values to establish MC as truth. GBBS was investigated for mono-energetic 1.0, 0.1, and 0.01 MeV electron and photon sources as well as 131I and 90Y radionuclides. We investigated convergence of GBBS by analyzing different meshes (), energy group structures () for each radionuclide component, angular quadrature orders (, and scattering order expansions (–); higher indices imply finer discretization. We compared GBBS to MC in (1) voxel-S-value geometry for soft tissue, lung, and bone, and (2) a source at the interface between combinations of lung, soft tissue, and bone.Excluding Auger and conversion electrons, MC agreed within ≈5% of published source voxel absorbed doses. For the finest discretization, most GBBS absorbed doses in the source voxel changed by less than 1% compared to the next finest discretization along each phase space variable indicating sufficient convergence. For the finest discretization, agreement with MC in the source voxel ranged from −3% to −20% with larger differences at lower energies (−3% for 1 MeV electron in lung to −20% for 0.01 MeV photon in bone); similar agreement was found for the interface geometries. Differences between GBBS and MC in the source voxel for 90Y and 131I were −6%.The GBBS ATTILA was benchmarked against MC in the nuclear medicine regime. GBBS can be a viable alternative to MC for voxel-level absorbed doses in nuclear medicine. However, reconciliation of the differences between GBBS and MC at lower energies requires further investigation of energy deposition cross-sections.

A NUMERICAL SCHEME FOR SPECIAL RELATIVISTIC RADIATION MAGNETOHYDRODYNAMICS BASED ON SOLVING THE TIME-DEPENDENT RADIATIVE TRANSFER EQUATION

ABSTRACT We develop a numerical scheme for solving the equations of fully special relativistic, radiation magnetohydrodynamics (MHDs), in which the frequency-integrated, time-dependent radiation transfer equation is solved to calculate the specific intensity. The radiation energy density, the radiation flux, and the radiation stress tensor are obtained by the angular quadrature of the intensity. In the present method, conservation of total mass, momentum, and energy of the radiation magnetofluids is guaranteed. We treat not only the isotropic scattering but also the Thomson scattering. The numerical method of MHDs is the same as that of our previous work. The advection terms are explicitly solved, and the source terms, which describe the gas–radiation interaction, are implicitly integrated. Our code is suitable for massive parallel computing. We present that our code shows reasonable results in some numerical tests for propagating radiation and radiation hydrodynamics. Particularly, the correct solution is given even in the optically very thin or moderately thin regimes, and the special relativistic effects are nicely reproduced.

Development and Evaluation of High-Fidelity Product and Evenly Spaced Angular Quadratures for Three-Dimensional Discrete Ordinates Calculations with Large Air Regions

Using two different methods, angular quadrature sets are developed with greater than about 1000 angles per octant to reduce ray effects in three-dimensional (3-D), discrete ordinates radiation transport calculations with large air or void regions. Quadrature sets from both methods are evaluated in two distinct 3-D models sensitive to quadrature details and are shown to behave reasonably well. The first method is a previously described method that is examined here in 3-D. The second method produces quadrature sets that have quadrature directions approximately evenly spaced over the entire surface of the unit sphere.

ACS Algorithm in Discrete Ordinates for Pressure Vessel Dosimetry

The Adaptive Collision Source (ACS) method can solve the Linear Boltzmann Equation (LBE) more efficiently by adaptation of the angular quadrature order. This is similar to, and essentially an extension of, the first collision source method. Previously, the ACS methodology has been implemented into the TITAN discrete ordinates code, and has shown speedups of 2-4 on a simple test problem, with very little loss of accuracy (within a provided adaptive tolerance). This work examines the use of the ACS method for a more realistic problem: pressure vessel dosimetry with the VENUS-2 MOX-fuelled reactor dosimetry benchmark. The ACS method proved to be able to obtain accurate results while being approximately twice as efficient as using a constant quadrature in a standard source iteration scheme.

Corner neutronic code

The neutronic code CORNER is based on the SN discrete ordinates method [1] and the PM scattering cross-section approximation. It is intended for high-precision deterministic neutronic calculations of fast-neutron reactors and can be used to solve two types of steady-state neutron transport and gamma quanta problems in a 3D hexagonal geometry: Keff problems (homogeneous) and source problems (inhomogeneous).The code is developed in the Fortran language and has a modular structure. Its key modules are a module for the preparation of neutron constants in the ANISN format; a geometric module containing a description of the core's loading map and fuel assembly types, including their axial meshing and material composition; a module for preparing angular quadrature sets; an input data module containing the parameters of the approximation used and the control parameters; a neutronic calculation module and a calculation data processing module.The Directional Theta-Weighted (DTW) difference scheme [2] has been built to approximate the spatial dependence. It has advantages over the DD (Diamond Difference) scheme broadly used in coarse-mesh problems.The energy dependence is represented by multigroup approximation. The angular variable is discretized by introduction of the angular quadrature set. Quadrature sets can be also defined by the user.An iterative solution process is used, including external iterations for the fission source and internal iterations for the scattering source. The paper presents the results of a cross-verification against the Monte Carlo MMK code [3] and on a model of the BN-800 reactor core.

Calculation of photon pulse height distribution using deterministic and Monte Carlo methods

Abstract Radiation transport techniques which are used in radiation detection systems comprise one of two categories namely probabilistic and deterministic. However, probabilistic methods are typically used in pulse height distribution simulation by recreating the behavior of each individual particle, the deterministic approach, which approximates the macroscopic behavior of particles by solution of Boltzmann transport equation, is being developed because of its potential advantages in computational efficiency for complex radiation detection problems. In current work linear transport equation is solved using two methods including collided components of the scalar flux algorithm which is applied by iterating on the scattering source and ANISN deterministic computer code. This approach is presented in one dimension with anisotropic scattering orders up to P8 and angular quadrature orders up to S16. Also, multi-group gamma cross-section library required for this numerical transport simulation is generated in a discrete appropriate form. Finally, photon pulse height distributions are indirectly calculated by deterministic methods that approvingly compare with those from Monte Carlo based codes namely MCNPX and FLUKA.

Advanced quadratures for three-dimensional discrete ordinate transport simulations: A comparative study

Three-dimensional radiation transport simulations play a key role in a number of problems found in nuclear science and engineering. When considering deterministic methods for solving the three-dimensional radiation transport equation, the Sn method of Carlson and Lee is one of the most commonly found techniques. The Sn method relies on a given set of “discrete ordinates” or discrete directions in which radiation streams. This set of discrete ordinates also must form a quadrature (numerical integration) so that the angular flux moments can be estimated and the scattering source calculated. In this paper we review a number of recently developed angular quadrature sets and compare their performance on a number of reactor physics problems.

The Source Equivalence Acceleration Method

Fine-group whole-core reactor analysis remains one of the long sought goals of the reactor physics community. Such a detailed analysis is typically too computationally expensive to be realized on anything except the largest of supercomputers. Recondensation using the Discrete Generalized Multigroup (DGM) method, though, offers a relatively cheap alternative to solving the fine group transport problem. DGM, however, suffered from inconsistencies when applied to high-order spatial methods. While an exact spatial recondensation method was developed and provided full spatial consistency with the fine group problem, this approach substantially increased memory requirements for realistic problems. The method described in this paper, called the Source Equivalence Acceleration Method (SEAM), forms a coarse-group problem which preserves the fine-group problem even when using higher order spatial methods. SEAM allows recondensation to converge to the fine-group solution with minimal memory requirements and little additional overhead. This method also provides for consistency when using different spatial methods and angular quadratures between the coarse group and fine group problems. SEAM was implemented in OpenMOC, a 2D MOC code developed at MIT, and its performance tested against Coarse Mesh Finite Difference (CMFD) acceleration on the C5G7 benchmark problem and on a 361 group version of the problem. For extremely expensive transport calculations, SEAM was able to outperform CMFD, resulting in speed-ups of 20–45 relative to the normal power iteration calculation.

LINTRAN v2.0: A linearised vector radiative transfer model for efficient simulation of satellite-born nadir-viewing reflection measurements of cloudy atmospheres

Radiance measurements of solar radiation that is backscattered by the Earth׳s atmosphere or surface contain information about the atmospheric composition and the state of the Earth׳s surface. Retrieving such information from satellite-based observations in nadir geometry employs a radiative transfer forward model. The forward model simulates the observed quantity, aiming to reproduce the observation.LINTRAN v2.0 is a linearised vector radiative transfer forward model, employing forward-adjoint theory, that is capable of modelling cloud contaminated satellite observations and their derivatives with respect to the state of the atmosphere and the Earth׳s surface in a numerically efficient manner. A significant gain in efficiency with respect to its predecessor (LINTRAN v1.0) is achieved through a mathematical framework that combines an approximate iterative solving method using the forward-adjoint perturbation theory with separation of the first N orders of scattering from the diffuse intensity vector field. Contributions to the observable up to order of scattering N are recursively solved in an analytical manner. Contributions from higher orders of scattering are subsequently solved in a numerical manner, assuming that the intensity field varies linearly with the vertical coordinate within an optically homogeneous model layer. This method is implemented in LINTRAN v2.0, choosing N=2, within the general framework of forward-adjoint perturbation theory.This new approach allows us to decrease the number of model layers and the degree of angular quadrature within the numerical solver by a factor of 10 and 1.4 respectively, compared to the previous model version, assuming a homogeneous atmosphere loaded with scattering Mie particles (size parameter χ≈35). In this homogeneous atmosphere, the reduced discretisation sampling in turn reduces the numerical effort associated with the numerical matrix solver by a factor of 42 relative to the previous model version, without a loss in model accuracy.

Sub-cell balanced nodal expansion methods using S4 eigenfunctions for multi-group SN transport problems in slab geometry

A highly accurate S4 eigenfunction-based nodal method has been developed to solve multi-group discrete ordinate neutral particle transport problems with a linearly anisotropic scattering in slab geometry. The new method solves the even-parity form of discrete ordinates transport equation with an arbitrary SN order angular quadrature using two sub-cell balance equations and the S4 eigenfunctions of within-group transport equation. The four eigenfunctions from S4 approximation have been chosen as basis functions for the spatial expansion of the angular flux in each mesh. The constant and cubic polynomial approximations are adopted for the scattering source terms from other energy groups and fission source. A nodal method using the conventional polynomial expansion and the sub-cell balances was also developed to be used for demonstrating the high accuracy of the new methods. Using the new methods, a multi-group eigenvalue problem has been solved as well as fixed source problems. The numerical test results of one-group problem show that the new method has third-order accuracy as mesh size is finely refined and it has much higher accuracies for large meshes than the diamond differencing method and the nodal method using sub-cell balances and polynomial expansion of angular flux. For multi-group problems including eigenvalue problem, it was demonstrated that the new method using the cubic polynomial approximation of the sources could produce very accurate solutions even with large mesh sizes.

Unique formulations in TITAN and PENTRAN for medical physics applications

Special algorithms for application to medical physical problems have been developed and implemented into the TITAN and PENTRAN 3-D parallel transport codes. TITAN is a 3-D parallel transport code with hybrid algorithms including 1) Sn and Characteristics, and 2) Sn with a fictitious quadrature set and ray-tracing. PENTRAN is a 3-D parallel transport code with adaptive differencing Sn formulation, full domain decomposition, and different angular quadrature types (including a characteristic option) with ordinate splitting. In this paper, we will discuss new algorithms developed in TITAN for image reconstruction in Single Photon Emission Computed Tomography (SPECT). It is demonstrated that TITAN projection images are in excellent agreement with the Monte Carlo MCNP5 predictions, while its computation time was up to 2787 times faster than the MCNP5 computation time for a 0.98° collimator acceptance angle. We will also discuss a novel algorithm developed in PENTRAN for dose calculations in heterogeneous, voxelized phantoms or other geometries. The new algorithm called EDK-Sn, or “Electron Dose Kernel-Discrete Ordinates (EDK-Sn)” yielded a speedup of ~8 over traditional highly parallel Monte Carlo calculation times, with <7% difference in dose among different organs (or smaller given stochastic uncertainties).

Revisiting the Fourier expansion of Mie scattering matrices in generalized spherical functions

Mie computations of the scattering properties of large particles are a time consuming step in the radiative transfer modeling of aerosol and clouds. Currently, there exist two methods based on the use of spherical functions for computing the Fourier moments of the phase matrix of a given spherical particle or particulate polydispersion: The first, developed over the years before being presented in a convenient form by Siewert [31], required an intermediate computation of the phase matrix over which numerical integration was performed to deliver the required Fourier components. The second, suggested by Domke [9], promised a direct computation of the Fourier moments using Wigner 3-j symbols.While the former was relatively easy to implement and is thus the most commonly used to date, its numerical implementation using an arbitrary user choice of angular quadrature (NAI-1) can produce inaccurate results. Numerical integration using quadrature points as recommended by de Rooij and van der Stap [5] (NAI-2) delivers accurate results with high computational efficiency.Domke's method enables a direct computation of the exact number of required Fourier components. However, the original manuscript contained several misprints, many of which were subsequently corrected by de Rooij and van der Stap [5]. Unfortunately, the main recurrence relationship used in Domke [9] remained uncorrected. In this paper, the corrected relationship is presented along with other minor corrections.de Rooij and van der Stap [5] had found the straightforward application of Domke's method viable only for size parameters smaller than ~120 due to issues involving computer storage. A means of implementing the corrected Domke formalism using precomputed tabulations of Wigner 3-j symbols (PCW) is presented here, making it more computationally economical and applicable over much broader particle size ranges. The accuracy of PCW is only limited by machine precision. For a single particle, NAI-2 is found to be faster than PCW for size parameters greater than about 228, whereas for polydispersions over a finite range of particle sizes, PCW is found to be at least 6–8 times faster for size parameters ranging from ~0 to beyond 900. PCW thus allows for significant reduction of the computational burden associated with Mie calculations for polydispersions.

Monte Carlo and Deterministic Neutronics Analyses of YALINA Thermal Facility and Comparison with Experimental Results

This study compares Monte Carlo and deterministic neutronics analyses of the zero-power YALINA Thermal subcritical assembly, which is located in Minsk, Belarus. The YALINA Thermal facility consists of a subcritical core that can be driven by either a californium neutron source or a deuterium-deuterium (D-D) neutron source. The californium neutron source is generated by the natural decay of 252Cf; the D-D neutron source is generated by a deuteron accelerator. The MCNPX, MONK, NJOY, DRAGON, PARTISN, and TORT computer programs have been used for calculating the neutron spectrum, the neutron flux, and the 3He(n,p) reaction rate set by californium and D-D neutron sources. These parameters have been computed in different experimental channels of the assembly for different fuel loading configurations. The MCNPX and MONK computer programs modeled the facility without any major approximation; the PARTISN and TORT computer simulations used 69 energy groups, S16 angular quadrature set, linear anisotropic scattering, and approximately 60 homogenized material zones. The results calculated by different computer programs are in good agreement; in addition, they match the 3He(n,p) reaction rate from experimental measurements obtained by californium and D-D neutron sources.

Analytical methods for computational modeling of fixed-source slab-geometry discrete ordinates transport problems: Response matrix and hybrid SN

Presented here are two analytical methods for computational modeling of fixed-source slab-geometry discrete ordinates (SN) transport problems for shielding and nuclear reactor fuel-moderator lattice calculations. For shielding calculations a response matrix method is described, which generates numerical solutions completely free from spatial truncation errors. Therefore, the offered response matrix method with the one-region block inversion (RBI) iterative scheme converges numerical results for the region-edge angular fluxes, that coincide with the numerical values generated from the SN analytical solution, apart from computational finite arithmetic considerations. As with lattice cell calculations in nuclear reactor physics, we describe an analytical direct method for hybrid SN calculations. The basic idea is to use higher order angular quadrature set in the highly absorbing fuel region (SNF) and lower order angular quadrature set in the diffusive moderator region (SNM), i.e., NF > NM. Special continuity conditions for the fuel exiting fluxes that constitute the incoming fluxes for the moderation region, and conversely, for the moderator exiting fluxes that constitute the incoming fluxes for the fuel region are applied in the present hybrid numerical scheme. These special continuity conditions are based on the equivalence of the multigroup SN and PN − 1 equations in slab geometry. A spatial reconstruction scheme is added to yield detailed profile of the solution within each homogenized region of the slab, starting from the coarse-mesh results. Test problems are given to illustrate the methods' accuracy.

A Comparison of Angular Discretization Strategies for Modeling Radiative Transfer in Pool Fire Simulations

Decoupled radiative heat transfer calculations of 30-cm-diameter heptane pool fire are carried out employing the P1, discrete ordinates (DO), and finite-volume (FV) methods. The DO calculations are performed employing three different quadrature schemes (level symmetric hybrid [LSH], spherical surface symmetrical equal dividing angular quadrature [SSD], and TN). By scaling the soot field in the fire, the performances of the various methods are evaluated at different optical thicknesses by employing highly angularly resolved FV calculations (512 directions) as a reference. Twenty-four angular directions in the DO and FV methods are seen to be adequate to accurately predict the distributions of the divergence of the radiative flux within the fire. The corresponding accuracies of the P1 model are seen to decrease with increase in optical thickness. In the regions away from the fire, the P1 model performs poorly in predicting the axial as well as the radial radiative fluxes, whereas the errors associated with employing coarse angular resolutions in the DO method are seen to increase with soot concentrations and distance from the pool centre. Among the DO and FV models, the accuracies of the axial and radial incident radiative flux predictions outside the fire were more dependent on the angular resolution of the calculations rather than on the particular quadrature set employed in the calculations.

A C5 Benchmark Problem with the Discrete Ordinates Radiation Transport Code Denovo

The C5 benchmark problem proposed by the Organisation for Economic Co-operation and Development/Nuclear Energy Agency was modeled to examine the capabilities of Denovo, a three-dimensional (3-D) parallel discrete ordinates (SN) radiation transport code, for problems with no spatial homogenization. Denovo uses state-of-the-art numerical methods to obtain accurate solutions to the Boltzmann transport equation. Problems were run in parallel on Jaguar, a high-performance supercomputer located at Oak Ridge National Laboratory. Both the two-dimensional (2-D) and 3-D configurations were analyzed, and the results were compared with the reference MCNP Monte Carlo calculations. For an additional comparison, SCALE/KENO-V.a Monte Carlo solutions were also included. In addition, a sensitivity analysis was performed for the optimal angular quadrature and mesh resolution for both the 2-D and 3-D infinite lattices of UO2 fuel pin cells. Denovo was verified with the C5 problem. The effective multiplication factors, pin powers, and assembly powers were found to be in good agreement with the reference MCNP and SCALE/KENO-V.a Monte Carlo calculations.