Stochastic methods for uncertainty quantification in radiation transport

The use of stochastic spectral expansions, specifically generalized polynomial chaos (gPC) and Karhunen-Loeve (KL) expansions, is investigated for uncertainty quantification in radiation transport. The gPC represents second-order random processes in terms of an expansion of orthogonal polynomials of random variables. The KL expansion is a Fourier-type expansion that represents a second-order random process with known covariance function in terms of a set of uncorrelated random variables and the eigenmodes of the covariance function. The flux and, in multiplying materials, the k-eigenvalue, which are the problem unknowns, are always expanded in a gPC expansion since their covariance functions are also unknown. This work assumes a single uncertain input—the total macroscopic cross section—although this does not represent a limitation of the approaches considered here. Two particular types of input parameter uncertainty are investigated: The cross section as a univariate Gaussian, log-normal, gamma or beta random variable, and the cross section as a spatially varying Gaussian or log-normal random process. In the first case, a gPC expansion in terms of a univariate random variable suffices, while in the second, a truncated KL expansion is first necessary followed by a gPC expansion in terms of multivariate random variables. Two solution methods are examined: The Stochastic Finite Element Method (SFEM) and the Stochastic Collocation Method (SCM). The SFEM entails taking Galerkin projections onto the orthogonal basis, which yields a system of fully-coupled equations for the PC coefficients of the flux and the k-eigenvalue. This system is linear when there is no multiplication and can be solved using Richardson iteration, employing a standard operator splitting such as block Gauss-Seidel or block Jacobi, or a Krylov iterative method, which can be preconditioned using these splittings. When multiplication is present, the SFEM system is non-linear and a Newton-Krylov method is employed to solve it. The SCM utilizes a suitable quadrature rule to compute the moments or PC coefficients of the flux and k-eigenvalue, and thus involves the solution of a system of independent deterministic transport equations. The accuracy and efficiency of the two methods are compared and contrasted. Both are shown to accurately compute the PC coefficients of the unknown, and numerical proof is provided that the two methods are in fact equivalent in certain cases. The PC coefficients are used to compute the moments and probability density functions of the unknowns, which are shown to be accurate by comparing with Monte Carlo results. An analytic diffusion analysis, corroborated by numerical results, reveals that the random transport equation is well approximated by a deterministic diffusion equation when the medium is diffusive with respect to the average cross section but without constraint on the amplitude of the random fluctuations. Our work shows that stochastic spectral expansions are a viable alternative to random sampling-based uncertainty quantification techniques since both provide a complete characterization of the distribution of the flux and the k-eigenvalue. Furthermore, it is demonstrated that, unlike perturbation methods, SFEM and SCM can handle large parameter uncertainty

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

Efficient fast Fourier transform-based solvers for computing the thermomechanical behavior of applied materials

The mechanical behavior of many applied materials arises from their microstructure. Thus, to aid the design, development and industrialization of new materials, robust computational homogenization methods are indispensable. The present thesis is devoted to investigating and developing FFT-based micromechanics solvers for efficiently computing the (thermo)mechanical response of nonlinear composite materials with complex microstructures

Coupling schemes and inexact Newton for multi-physics and coupled optimization problems

This work targets mathematical solutions and software for complex numerical simulation and optimization problems. Characteristics are the combination of different models and software modules and the need for massively parallel execution on supercomputers. We consider two different types of multi-component problems in Part I and Part II of the thesis: (i) Surface coupled fluid- structure interactions and (ii) analysis of medical MR imaging data of brain tumor patients. In (i), we establish highly accurate simulations by combining different aspects such as fluid flow and arterial wall deformation in hemodynamics simulations or fluid flow, heat transfer and mechanical stresses in cooling systems. For (ii), we focus on (a) facilitating the transfer of information such as functional brain regions from a statistical healthy atlas brain to the individual patient brain (which is topologically different due to the tumor), and (b) to allow for patient specific tumor progression simulations based on the estimation of biophysical parameters via inverse tumor growth simulation (given a single snapshot in time, only). Applications and specific characteristics of both problems are very distinct, yet both are hallmarked by strong inter-component relations and result in formidable, very large, coupled systems of partial differential equations. Part I targets robust and efficient quasi-Newton methods for black-box surface-coupling of parti- tioned fluid-structure interaction simulations. The partitioned approach allows for great flexibility and exchangeable of sub-components. However, breaking up multi-physics into single components requires advanced coupling strategies to ensure correct inter-component relations and effectively tackle instabilities. Due to the black-box paradigm, solver internals are hidden and information exchange is reduced to input/output relations. We develop advanced quasi-Newton methods that effectively establish the equation coupling of two (or more) solvers based on solving a non-linear fixed-point equation at the interface. Established state of the art methods fall short by either requiring costly tuning of problem dependent parameters, or becoming infeasible for large scale problems. In developing parameter-free, linear-complexity alternatives, we lift the robustness and parallel scalability of quasi-Newton methods for partitioned surface-coupled multi-physics simulations to a new level. The developed methods are implemented in the parallel, general purpose coupling tool preCICE. Part II targets MR image analysis of glioblastoma multiforme pathologies and patient specific simulation of brain tumor progression. We apply a joint medical image registration and biophysical inversion strategy, targeting at facilitating diagnosis, aiding and supporting surgical planning, and improving the efficacy of brain tumor therapy. We propose two problem formulations and decompose the resulting large-scale, highly non-linear and non-convex PDE-constrained optimization problem into two tightly coupled problems: inverse tumor simulation and medical image registration. We deduce a novel, modular Picard iteration-type solution strategy. We are the first to successfully solve the inverse tumor-growth problem based on a single patient snapshot with a gradient-based approach. We present the joint inversion framework SIBIA, which scales to very high image resolutions and parallel execution on tens of thousands of cores. We apply our methodology to synthetic and actual clinical data sets and achieve excellent normal-to-abnormal registration quality and present a proof of concept for a very promising strategy to obtain clinically relevant biophysical information. Advanced inexact-Newton methods are an essential tool for both parts. We connect the two parts by pointing out commonalities and differences of variants used in the two communities in unified notation

Efficient fast Fourier transform-based solvers for computing the thermomechanical behavior of applied materials

The mechanical behavior of many applied materials arises from their microstructure. Thus, to aid the design, development and industrialization of new materials, robust computational homogenization methods are indispensable. The present thesis is devoted to investigating and developing FFT-based micromechanics solvers for efficiently computing the (thermo)mechanical response of nonlinear composite materials with complex microstructures

Mathematical and Algorithmic Aspects of Atmosphere-Ocean Data Assimilation

The nomenclature “data assimilation” arises from applications in the geosciences where complex mathematical models are interfaced with observational data in order to improve model forecasts. Mathematically, data assimilation is closely related to filtering and smoothing on the one hand and inverse problems and statistical inference on the other. Key challenges of data assimilation arise from the high-dimensionality of the underlying models, combined with systematic spatio-temporal model errors, pure model uncertainty quantifications and relatively sparse observation networks. Advances in the field of data assimilation will require combination of a broad range of mathematical techniques from differential equations, statistics, probability, scientific computing and mathematical modelling, together with insights from practitioners in the field. The workshop brought together a collection of scientists representing this broad spectrum of research strands

Tomographic Inverse Problems: Theory and Applications

This was the tenth Oberwolfach conference on the mathematics of tomography. The field rests on the interplay between the theoretical and applied; practical questions lead to new mathematics and pure mathematics motivates new algorithms. This workshop encompassed classical areas such as X-ray computed tomography (CT) as well as new modalities and applications such as dynamic imaging, Compton scattering tomography, hybrid imaging, optical tomography or multi-energy CT and addressed inter alia the use of methods from machine learning

Hybrid Multiscale Methods with Applications to Semiconductors, Porous Media, and Materials Science

In this work we consider two multiscale applications with tremendous computational complexity at the lower scale. First, we examine a model for charge transport in semicon- ductor structures with heterojunction interfaces. Due to the complex physical phenomena at the interface, the model at the design scale is unable to adequately capture the behavior of the structure in the interface region. Simultaneously it is computationally intractable to simulate the full heterostructure on the scale required near the interface. Second, we con- sider the problem of the simulation of fluid flow in a dynamically evolving porous medium. The evolution of the medium strongly couples the porescale flow solutions and the macro scale model, requiring a novel approach to communicate the porescale evolution to the macroscale without resorting to the intractable simulation of the fluid flow problem di- rectly on the porescale geometry. We formulate novel methods for these two applications in the multiscale framework. For the semiconductor problem we present iterative sub- structuring domain decomposition methods that decouple the interface computation from the macroscale model. For the fluid flow problem we develop a reduced order three-scale fluid flow model based on a spatial decomposition of the porescale geometry and the offline approximation of a stochastic process describing macroscale permeability paramaterized by the volume fraction of the evolved geometry