Improved Fully-Implicit Spherical Harmonics Methods for First and Second Order Forms of the Transport Equation Using Galerkin Finite Element

In this dissertation, we focus on solving the linear Boltzmann equation -- or transport equation -- using spherical harmonics (PN) expansions with fully-implicit time-integration schemes and Galerkin Finite Element spatial discretizations within the Multiphysics Object Oriented Simulation Environment (MOOSE) framework. The presentation is composed of two main ensembles. On one hand, we study the first-order form of the transport equation in the context of Thermal Radiation Transport (TRT). This nonlinear application physically necessitates to maintain a positive material temperature while the PN approximation tends to create oscillations and negativity in the solution. To mitigate these flaws, we provide a fully-implicit implementation of the Filtered PN (FPN) method and investigate local filtering strategies. After analyzing its effect on the conditioning of the system and showing that it improves the convergence properties of the iterative solver, we numerically investigate the error estimates derived in the linear setting and observe that they hold in the non-linear case. Then, we illustrate the benefits of the method on a standard test problem and compare it with implicit Monte Carlo (IMC) simulations. On the other hand, we focus on second-order forms of the transport equation for neutronics applications. We mostly consider the Self-Adjoint Angular Flux (SAAF) and Least-Squares (LS) formulations, the former being globally conservative but void incompatible and the latter having -- in all generality -- the opposite properties. We study the relationship between these two methods based on the weakly-imposed LS boundary conditions. Equivalences between various parity-based PN methods are also established, in particular showing that second-order filters are not an appropriate fix to retrieve void compatibility. The importance of global conservation is highlighted on a heterogeneous multigroup k-eigenvalue test problem. Based on these considerations, we propose a new method that is both globally conservative and compatible with voids. The main idea is to solve the LS form in the void regions and the SAAF form elsewhere. For the LS form to be conservative in void, a non-symmetric fix is required, yielding the Conservative LS (CLS) formulation. From there, a hybrid SAAF-- CLS method can be derived, having the desired properties. We also show how to extend it to near-void regions and time-dependent problems. While such a second-order form already existed for discrete-ordinates (SN) discretizations (Wang et al. 2014), we believe that this method is the first of its kind, being well-suited to both SN and PN discretizations

Multilevel Schwarz methods for multigroup radiation transport problems

The development of advanced discretization methods for the radiation transport equation is of fundamental importance, since the numerical effort of modeling increasingly complex multidimensional problems with increasing accuracy is extremely challenging. Different expressions of this equation arise in several science fields, from nuclear fission and fusion to astrophysics, climatology and combustion. Mathematically, the radiation intensity is usually a rapidly changing function, causing a considerable loss in accuracy for many discretization methods. Depending on the coefficient ranges, the equation behaves like totally different equation types, making it very difficult to find a discretization method that is efficient in all regimes. Computationally, the huge amount of unknowns involved demands not only extremely powerful computers, but also efficient numerical methods and optimized implementations. Today, solvers covering all the coefficient ranges and still being robust in the diffusion dominated case are very scarce. In the last 20 years, Discontinous Galerkin (DG) methods have been studied for the monoenergetic problem, unsuccessfully, due to lack of stability for diffusion-dominated cases. Recently, new mathematical developments have fully explained the instability and provided a remedy by using a numerical flux depending on the scattering cross section and the mesh size. The new formulation has proven to be stable and allows the application of multigrid, matrix-free methods, reducing the memory needed for such an amount of unknowns. We use these numerical methods to address the solution of a energy dependent problem with a multigroup approach. We study the diffusion approximation to the transport problem, obtaining convergence proofs for the symmetric scattering case and advances in the nonsymmetric case, using field of values analysis. For the full transport case, we discretize by means of an asymptotic preserving, weakly penalized discontinuous Galerkin method that we solve with a multigrid preconditioned GMRES solver, using nonoverlapping Schwarz smoothers for the energy and direction dependent radiative transfer problem. To address the local thermodynamic equilibrium (LTE) constraint, we use a nonlinear additive Schwarz method to precondition the Newton solver. By solving full local radiative transfer problems for each grid cell, performed in parallel on a matrix-free implementation, we achieve a method capable to address large scale calculations arising from applications such as astrophysics, atmospheric radiation calculations and nuclear applications. To the best of our knowledge, this is the first time this preconditioner combination has been used in LTE radiation transport and in several tests we show the robustness of the approach for different mesh sizes, cross sections, energy distributions and anisotropic regimes, both in the linear and nonlinear cases

Adaptive Mesh Refinement Solution Techniques for the Multigroup SN Transport Equation Using a Higher-Order Discontinuous Finite Element Method

In this dissertation, we develop Adaptive Mesh Refinement (AMR) techniques for the steady-state multigroup SN neutron transport equation using a higher-order Discontinuous Galerkin Finite Element Method (DGFEM). We propose two error estimations, a projection-based estimator and a jump-based indicator, both of which are shown to reliably drive the spatial discretization error down using h-type AMR. Algorithms to treat the mesh irregularity resulting from the local refinement are implemented in a matrix-free fashion. The DGFEM spatial discretization scheme employed in this research allows the easy use of adapted meshes and can, therefore, follow the physics tightly by generating group-dependent adapted meshes. Indeed, the spatial discretization error is controlled with AMR for the entire multigroup SNtransport simulation, resulting in group-dependent AMR meshes. The computing efforts, both in memory and CPU-time, are significantly reduced. While the convergence rates obtained using uniform mesh refinement are limited by the singularity index of transport solution (3/2 when the solution is continuous, 1/2 when it is discontinuous), the convergence rates achieved with mesh adaptivity are superior. The accuracy in the AMR solution reaches a level where the solution angular error (or ray effects) are highlighted by the mesh adaptivity process. The superiority of higherorder calculations based on a matrix-free scheme is verified on modern computing architectures. A stable symmetric positive definite Diffusion Synthetic Acceleration (DSA) scheme is devised for the DGFEM-discretized transport equation using a variational argument. The Modified Interior Penalty (MIP) diffusion form used to accelerate the SN transport solves has been obtained directly from the DGFEM variational form of the SN equations. This MIP form is stable and compatible with AMR meshes. Because this MIP form is based on a DGFEM formulation as well, it avoids the costly continuity requirements of continuous finite elements. It has been used as a preconditioner for both the standard source iteration and the GMRes solution technique employed when solving the transport equation. The variational argument used in devising transport acceleration schemes is a powerful tool for obtaining transportconforming diffusion schemes. xuthus, a 2-D AMR transport code implementing these findings, has been developed for unstructured triangular meshes

Crash: A Block-Adaptive-Mesh Code for Radiative Shock Hydrodynamics - Implementation and Verification

We describe the CRASH (Center for Radiative Shock Hydrodynamics) code, a block adaptive mesh code for multi-material radiation hydrodynamics. The implementation solves the radiation diffusion model with the gray or multigroup method and uses a flux limited diffusion approximation to recover the free-streaming limit. The electrons and ions are allowed to have different temperatures and we include a flux limited electron heat conduction. The radiation hydrodynamic equations are solved in the Eulerian frame by means of a conservative finite volume discretization in either one, two, or three-dimensional slab geometry or in two-dimensional cylindrical symmetry. An operator split method is used to solve these equations in three substeps: (1) solve the hydrodynamic equations with shock-capturing schemes, (2) a linear advection of the radiation in frequency-logarithm space, and (3) an implicit solve of the stiff radiation diffusion, heat conduction, and energy exchange. We present a suite of verification test problems to demonstrate the accuracy and performance of the algorithms. The CRASH code is an extension of the Block-Adaptive Tree Solarwind Roe Upwind Scheme (BATS-R-US) code with this new radiation transfer and heat conduction library and equation-of-state and multigroup opacity solvers. Both CRASH and BATS-R-US are part of the publicly available Space Weather Modeling Framework (SWMF).Comment: 51 pages, 19 figures; submitted to Astrophysical Journa

An adaptive Cartesian embedded boundary approach for fluid simulations of two- and three-dimensional low temperature plasma filaments in complex geometries

We review a scalable two- and three-dimensional computer code for low-temperature plasma simulations in multi-material complex geometries. Our approach is based on embedded boundary (EB) finite volume discretizations of the minimal fluid-plasma model on adaptive Cartesian grids, extended to also account for charging of insulating surfaces. We discuss the spatial and temporal discretization methods, and show that the resulting overall method is second order convergent, monotone, and conservative (for smooth solutions). Weak scalability with parallel efficiencies over 70\% are demonstrated up to 8192 cores and more than one billion cells. We then demonstrate the use of adaptive mesh refinement in multiple two- and three-dimensional simulation examples at modest cores counts. The examples include two-dimensional simulations of surface streamers along insulators with surface roughness; fully three-dimensional simulations of filaments in experimentally realizable pin-plane geometries, and three-dimensional simulations of positive plasma discharges in multi-material complex geometries. The largest computational example uses up to $800$ million mesh cells with billions of unknowns on $4096$ computing cores. Our use of computer-aided design (CAD) and constructive solid geometry (CSG) combined with capabilities for parallel computing offers possibilities for performing three-dimensional transient plasma-fluid simulations, also in multi-material complex geometries at moderate pressures and comparatively large scale.Comment: 40 pages, 21 figure

Higher-Order DGFEM Transport Calculations on Polytope Meshes for Massively-Parallel Architectures

In this dissertation, we develop improvements to the discrete ordinates (S_N) neutron transport equation using a Discontinuous Galerkin Finite Element Method (DGFEM) spatial discretization on arbitrary polytope (polygonal and polyhedral) grids compatible for massively-parallel computer architectures. Polytope meshes are attractive for multiple reasons, including their use in other physics communities and their ease in handling local mesh refinement strategies. In this work, we focus on two topical areas of research. First, we discuss higher-order basis functions compatible to solve the DGFEM S_N transport equation on arbitrary polygonal meshes. Second, we assess Diffusion Synthetic Acceleration (DSA) schemes compatible with polytope grids for massively-parallel transport problems. We first utilize basis functions compatible with arbitrary polygonal grids for the DGFEM transport equation. We analyze four different basis functions that have linear completeness on polygons: the Wachspress rational functions, the PWL functions, the mean value coordinates, and the maximum entropy coordinates. We then describe the procedure to extend these polygonal linear basis functions into the quadratic serendipity space of functions. These quadratic basis functions can exactly interpolate monomial functions up to order 2. Both the linear and quadratic sets of basis functions preserve transport solutions in the thick diffusion limit. Maximum convergence rates of 2 and 3 are observed for regular transport solutions for the linear and quadratic basis functions, respectively. For problems that are limited by the regularity of the transport solution, convergence rates of 3/2 (when the solution is continuous) and 1/2 (when the solution is discontinuous) are observed. Spatial Adaptive Mesh Refinement (AMR) achieved superior convergence rates than uniform refinement, even for problems bounded by the solution regularity. We demonstrated accuracy in the AMR solutions by allowing them to reach a level where the ray effects of the angular discretization are realized. Next, we analyzed DSA schemes to accelerate both the within-group iterations as well as the thermal upscattering iterations for multigroup transport problems. Accelerating the thermal upscattering iterations is important for materials (e.g., graphite) with significant thermal energy scattering and minimal absorption. All of the acceleration schemes analyzed use a DGFEM discretization of the diffusion equation that is compatible with arbitrary polytope meshes: the Modified Interior Penalty Method (MIP). MIP uses the same DGFEM discretization as the transport equation. The MIP form is Symmetric Positive De_nite (SPD) and e_ciently solved with Preconditioned Conjugate Gradient (PCG) with Algebraic MultiGrid (AMG) preconditioning. The analysis from previous work was extended to show MIP's stability and robustness for accelerating 3D transport problems. MIP DSA preconditioning was implemented in the Parallel Deterministic Transport (PDT) code at Texas A&M University and linked with the HYPRE suite of linear solvers. Good scalability was numerically verified out to around 131K processors. The fraction of time spent performing DSA operations was small for problems with sufficient work performed in the transport sweep (O(10^3) angular directions). Finally, we have developed a novel methodology to accelerate transport problems dominated by thermal neutron upscattering. Compared to historical upscatter acceleration methods, our method is parallelizable and amenable to massively parallel transport calculations. Speedup factors of about 3-4 were observed with our new method

Perturbed, Entropy-Based Closure for Radiative Transfer

We derive a hierarchy of closures based on perturbations of well-known entropy-based closures; we therefore refer to them as perturbed entropy-based models. Our derivation reveals final equations containing an additional convective and diffusive term which are added to the flux term of the standard closure. We present numerical simulations for the simplest member of the hierarchy, the perturbed M1 or PM1 model, in one spatial dimension. Simulations are performed using a Runge-Kutta discontinuous Galerkin method with special limiters that guarantee the realizability of the moment variables and the positivity of the material temperature. Improvements to the standard M1 model are observed in cases where unphysical shocks develop in the M1 model.Comment: 35 pages, 8 figure