In preparation for, and for decision support during, CBRN (chemical, biological, radiological and nuclear) emergencies it is essential to know how such an event would turn out, so that one can prepare a possible evacuation. Afterwards it might be good to know how to backtrack and see what caused the emergency, and in the case of e.g. a gas leak, where did it begin? The Swedish Defence Research Agency (FOI) develops models for such scenarios. In this thesis FOI's model, "The Nordic Urban Dispersion model" (NUD), has been studied. The system of equations set up by this model was originally solved using Intel's PARDISO solver, which is a direct solver. An evaluation on how an iterative multigrid method would work to solve the system has been done in this thesis. The wind model is a mass-consistent model which sets up a diagnostic initial wind field. The final wind field is later minimized under the constraint of the continuity equation. The minimization problem is solved using Lagrange multipliers and the system turns into a Poisson-like problem. The iterative algebraic multigrid solver (AMG) which has been evaluated had difficulties solving the problem of an asymmetric system matrix generated by NUD. The AMG solver was then tried on a symmetric discrete Poisson problem instead, and the solution turns out to be the same as for the PARDISO solver. A comparison was made between the AMG and PARDISO solver, and for the discrete Poisson case the AMG solver turned out on top for both larger system size and less computational time. To try out the solvers for the original NUD case a modification of the boundary conditions was made to make the system matrix symmetric. This modification turns the problem into a mathematical problem rather than a physical one, as the wind fields generated are not physically correct. For this modified case both the solvers get the same solution in essentially the same computational time. A method of how to in the future solve the original (asymmetric) problem, by modifying the discretization of the boundary conditions, has been discussed.

As finite element discretizations ever grow in size to address real-world problems, there is an increasing need for fast algorithms. Nowadays there are many GPU/CPU parallel approaches to solve such problems. Multigrid methods can be used to solve large-scale problems, or even better they can be used to precondition the conjugate gradient method, yielding better results in general. Capabilities of multigrid algorithms rely on the effectiveness of the inter-grid transfer operators. In this thesis we focus on the aggregation approach, discussing how different aggregation strategies affect the convergence rate. Based on these discussions, we propose an alternative parallel aggregation algorithm to improve convergence. We also provide numerous experimental results that compare different aggregation approaches, multigrid methods, and conjugate gradient iteration counts, showing that our proposed algorithm performs better in serial and parallel. / Modeling real-world problems incurs a high computational cost because these mathematical models involve large-scale data manipulation. Thus we need fast and efficient algorithms. Nowadays there are many high-performance approaches for these problems. One such method is called the Multigrid algorithm. This approach models a physical domain using a hierarchy of grids, and so the effectiveness of these approaches relies on how well data can be transferred from grid to grid. In this thesis, we focus on the aggregation approach, which clusters a grid’s vertices according to its connections. We also provide an alternative parallel aggregation algorithm to give a faster solution. We show numerous experimental results that compare different aggregation approaches and multigrid methods, showing that our proposed algorithm performs better in serial and parallel than other popular implementations.

The solution of large and sparse eigenvalue problems arising from the discretization of the diffusion equation is considered. The multigroup diffusion equation is discretized by means of the Nodal expansion Method (NEM) [9, 10]. A new formulation of the higher order NEM variants revealing the true nature of the problem, that is, a generalized eigenvalue problem, is proposed. These generalized eigenvalue problems are solved using the Jacobi-Davidson (JD) method [26]. The most expensive part of the method consists of solving a linear system referred to as correction equation. It is solved using Krylov subspace methods in combination with aggregation-based Algebraic Multigrid (AMG) techniques. In that context, a particular aggregation technique used in combination with classical smoothers, referred to as oblique geometric coarsening, has been derived. Its particularity is that it aggregates unknowns that are not coupled, which has never been done to our knowledge. A modular code, combining JD with an AMG preconditioner, has been developed. The code comes with many options, that have been tested. In particular, the instability of the Rayleigh-Ritz [33] acceleration procedure in the non-symmetric case has been underlined. Our code has also been compared to an industrial code extracted from ARTEMIS.

The majority of this thesis is concerned with the development of efficient and robust numerical methods based on adaptive algebraic multigrid to compute the stationary distribution of Markov chains. It is shown that classical algebraic multigrid techniques can be applied in an exact interpolation scheme framework to compute the stationary distribution of irreducible, homogeneous Markov chains. A quantitative analysis shows that algebraically smooth multiplicative error is locally constant along strong connections in a scaled system operator, which suggests that classical algebraic multigrid coarsening and interpolation can be applied to the class of nonsymmetric irreducible singular M-matrices with zero column sums. Acceleration schemes based on fine-level iterant recombination, and over-correction of the coarse-grid correction are developed to improve the rate of convergence and scalability of simple adaptive aggregation multigrid methods for Markov chains. Numerical tests over a wide range of challenging nonsymmetric test problems demonstrate the effectiveness of the proposed multilevel method and the acceleration schemes. This thesis also investigates the application of adaptive algebraic multigrid techniques for computing the canonical decomposition of higher-order tensors. The canonical decomposition is formulated as a least squares optimization problem, for which local minimizers are computed by solving the first-order optimality equations. The proposed multilevel method consists of two phases: an adaptive setup phase that uses a multiplicative correction scheme in conjunction with bootstrap algebraic multigrid interpolation to build the necessary operators on each level, and a solve phase that uses additive correction cycles based on the full approximation scheme to efficiently obtain an accurate solution. The alternating least squares method, which is a standard one-level iterative method for computing the canonical decomposition, is used as the relaxation scheme. Numerical tests show that for certain test problems arising from the discretization of high-dimensional partial differential equations on regular lattices the proposed multilevel method significantly outperforms the standard alternating least squares method when a high level of accuracy is required.

The majority of this thesis is concerned with the development of efficient and robust numerical methods based on adaptive algebraic multigrid to compute the stationary distribution of Markov chains. It is shown that classical algebraic multigrid techniques can be applied in an exact interpolation scheme framework to compute the stationary distribution of irreducible, homogeneous Markov chains. A quantitative analysis shows that algebraically smooth multiplicative error is locally constant along strong connections in a scaled system operator, which suggests that classical algebraic multigrid coarsening and interpolation can be applied to the class of nonsymmetric irreducible singular M-matrices with zero column sums. Acceleration schemes based on fine-level iterant recombination, and over-correction of the coarse-grid correction are developed to improve the rate of convergence and scalability of simple adaptive aggregation multigrid methods for Markov chains. Numerical tests over a wide range of challenging nonsymmetric test problems demonstrate the effectiveness of the proposed multilevel method and the acceleration schemes. This thesis also investigates the application of adaptive algebraic multigrid techniques for computing the canonical decomposition of higher-order tensors. The canonical decomposition is formulated as a least squares optimization problem, for which local minimizers are computed by solving the first-order optimality equations. The proposed multilevel method consists of two phases: an adaptive setup phase that uses a multiplicative correction scheme in conjunction with bootstrap algebraic multigrid interpolation to build the necessary operators on each level, and a solve phase that uses additive correction cycles based on the full approximation scheme to efficiently obtain an accurate solution. The alternating least squares method, which is a standard one-level iterative method for computing the canonical decomposition, is used as the relaxation scheme. Numerical tests show that for certain test problems arising from the discretization of high-dimensional partial differential equations on regular lattices the proposed multilevel method significantly outperforms the standard alternating least squares method when a high level of accuracy is required.

The development of adaptive mesh refinement capabilities in the field of computational fluid dynamics is an essential tool for enabling the simulation of larger and more complex physical problems. While such techniques have been known for a long time, most simulations do not make use of them because of the lack of a robust implementation. In this work, we present recent progresses that have been made to develop adaptive mesh refinement features in Nek5000, a code based on the spectral element method. These developments are driven by the algorithmic challenges posed by future exascale supercomputers. First, we perform the study of the strong scaling of Nek5000 on three petascale machines in order to assess the scalability of the code and identify the current bottlenecks. It is found that strong scaling limit ranges between 5, 000 and 220, 000 degrees of freedom per core depending on the machine and the case. The need for synchronized and low latency communication for efficient computational fluid dynamics simulation is also confirmed. Additionally, we present how Hypre, a library for linear algebra, is used to develop a new and efficient code for performing the setup step required prior to the use of an algebraic multigrid solver for preconditioning the pressure equation in Nek5000. Finally, the main objective of this work is to develop new methods for estimating the error on a numerical solution of the Navier–Stokes equations via the resolution of an adjoint problem. These new estimators are compared to existing ones, which are based on the decay of the spectral coefficients. Then, the estimators are combined with newly implemented capabilities in Nek5000 for automatic grid refinement and adaptive mesh adaptation is carried out. The applications considered so far are steady and two-dimensional, namely the lid-driven cavity at Re = 7, 500 and the flow past a cylinder at Re = 40. The use of adaptive mesh refinement techniques makes mesh generation easier and it is shown that a similar accuracy as with a static mesh can be reached with a significant reduction in the number of degrees of freedom. / <p>QC 20171114</p>

The implementation of a hybrid spectral/finite-element discretization on the unsteady, incompressible, Navier-Stokes equations with a semi-implicit time-stepping method, an explicit treatment of the advective terms, and an implicit treatment of the pressure and viscous terms leads to an algorithm capable of calculating 3D flows over complex 2D geometries. This also results in multiple Fourier mode linear systems which must be solved at every timestep, which naturally leads to two parallelization approaches: Fourier space partitioning, where each processor individually and simultaneously solves a linear system, and physical space partitioning, where all processors collectively solve each linear system, sequentially advancing through Fourier modes. These two parallelization approaches are compared based upon computational cost using multiple solvers: direct sparse LU, smoothed aggregation AMG, and single-level ILUT preconditioned GMRES; and on two supercomputers of different memory architecture(distributed and shared memory). This study revealed Fourier space partitioning outperforms physical space partitioning in all problems analyzed, and scales more efficiently as well. These differences were more dramatic on the distributed memory platform than the shared memory platform. Another study compares the previously mentioned solvers along with one additional solver, pointwise AMG, in Fourier space partitioning without parallelization to better understand computational scaling for problems with large meshes. It was found that the direct sparse LU solver performed well in terms of computational time, scaled linearly, but had very high memory usage which scaled in a super-linear manner. The single-level ILUT preconditioned GMRES solver required the least amount of memory, which also scaled linearly, but only had acceptable performance in terms of computational time for coarse meshes. Both AMG methods scaled linearly in both memory usage and time, and were comparable to the direct sparse LU solver in terms of computational time. The results of these studies are particularly useful for implementation of this algorithm on challenging and complex flows, especially direct numerical and large-eddy simulations. Reducing computational cost allows the analysis and understanding of more flows of practical interest.

Este trabalho consiste na análise do comportamento do método WAMG (Wavelet-Based Algebraic Multigrid), método numérico de resolução de sistemas de equações lineares desenvolvido no LMAG-Laboratório de Eletromagnetismo Aplicado, com relação a diversos suavizadores. O fato dos vetores que compõem os operadores matriciais Pronlongamento e Restrição do método WAMG serem ortonormais viabiliza uma série de análises teóricas e de dados experimentais, permitindo visualizar características não permitidas nos outros métodos Multigrid (MG), englobando o Multigrid Geométrico (GMG) e o Multigrid Algébrico (AMG). O método WAMG V-Cycle com Filtro Haar é testado em uma variedade de sistemas de equações lineares variando o suavizador, o coeficiente de relaxação nos suavizadores Damped Jacobi e Sobre Relaxação Sucessiva (SOR), e a configuração de pré e pós-suavização. Entre os suavizadores testados, estão os métodos iterativos estacionários Damped Jacobi, SOR, Esparsa Aproximada a Inversa tipo Diagonal (SPAI-0) e métodos propostos com a característica de suavização para-otimizada. A título de comparação, métodos iterativos não estacionários são testados também como suavizadores como Gradientes Conjugados, Gradientes Bi-Conjugados e ICCG. Os resultados dos testes são apresentados e comentados. / This work is comprised of WAMG (Wavelet-Based Algebraic Multigrid) method behavioral analysis based on variety of smoothers, numerical method based on linear equation systems resolution developed at LMAG (Applied Electromagnetism Laboratory). Based on the fact that the vectors represented by WAMG Prolongation and Restriction matrix operators are orthonormals allows the use of a variety of theoretical and practical analysis, and therefore gain visibility of characteristics not feasible through others Multigrid (MG) methods, such as Geometric Multigrid (GMG) and Algebraic Multigrid (AMG). WAMG V-Cycle method with Haar Filter is tested under a variety of linear equation systems, by varying smoothers, relaxation coefficient at Damped Jacobi and Successive Over Relaxation (SOR) smoothers, and pre and post smoothers configurations. The tested smoothers are stationary iterative methods such as Damped Jacobi, SOR, Diagonal type-Sparse Approximate Inverse (SPAI-0) and suggested ones with optimized smoothing characteristic. For comparison purposes, the Conjugate Gradients, Bi-Conjugate Gradient and ICCG non-stationary iterative methods are also tested as smoothers. The testing results are formally presented and commented.

Este trabalho consiste na análise do comportamento do método WAMG (Wavelet-Based Algebraic Multigrid), método numérico de resolução de sistemas de equações lineares desenvolvido no LMAG-Laboratório de Eletromagnetismo Aplicado, com relação a diversos suavizadores. O fato dos vetores que compõem os operadores matriciais Pronlongamento e Restrição do método WAMG serem ortonormais viabiliza uma série de análises teóricas e de dados experimentais, permitindo visualizar características não permitidas nos outros métodos Multigrid (MG), englobando o Multigrid Geométrico (GMG) e o Multigrid Algébrico (AMG). O método WAMG V-Cycle com Filtro Haar é testado em uma variedade de sistemas de equações lineares variando o suavizador, o coeficiente de relaxação nos suavizadores Damped Jacobi e Sobre Relaxação Sucessiva (SOR), e a configuração de pré e pós-suavização. Entre os suavizadores testados, estão os métodos iterativos estacionários Damped Jacobi, SOR, Esparsa Aproximada a Inversa tipo Diagonal (SPAI-0) e métodos propostos com a característica de suavização para-otimizada. A título de comparação, métodos iterativos não estacionários são testados também como suavizadores como Gradientes Conjugados, Gradientes Bi-Conjugados e ICCG. Os resultados dos testes são apresentados e comentados. / This work is comprised of WAMG (Wavelet-Based Algebraic Multigrid) method behavioral analysis based on variety of smoothers, numerical method based on linear equation systems resolution developed at LMAG (Applied Electromagnetism Laboratory). Based on the fact that the vectors represented by WAMG Prolongation and Restriction matrix operators are orthonormals allows the use of a variety of theoretical and practical analysis, and therefore gain visibility of characteristics not feasible through others Multigrid (MG) methods, such as Geometric Multigrid (GMG) and Algebraic Multigrid (AMG). WAMG V-Cycle method with Haar Filter is tested under a variety of linear equation systems, by varying smoothers, relaxation coefficient at Damped Jacobi and Successive Over Relaxation (SOR) smoothers, and pre and post smoothers configurations. The tested smoothers are stationary iterative methods such as Damped Jacobi, SOR, Diagonal type-Sparse Approximate Inverse (SPAI-0) and suggested ones with optimized smoothing characteristic. For comparison purposes, the Conjugate Gradients, Bi-Conjugate Gradient and ICCG non-stationary iterative methods are also tested as smoothers. The testing results are formally presented and commented.

Projecting the ice sheets' contribution to sea-level rise is difficult because of the complexity of accurately modeling ice sheet dynamics for the full polar ice sheets, because of the uncertainty in key, unobservable parameters governing those dynamics, and because quantifying the uncertainty in projections is necessary when determining the confidence to place in them. This work presents the formulation and solution of the Bayesian inverse problem of inferring, from observations, a probability distribution for the basal sliding parameter field beneath the Antarctic ice sheet. The basal sliding parameter is used within a high-fidelity nonlinear Stokes model of ice sheet dynamics. This model maps the parameters "forward" onto a velocity field that is compared against observations. Due to the continental-scale of the model, both the parameter field and the state variables of the forward problem have a large number of degrees of freedom: we consider discretizations in which the parameter has more than 1 million degrees of freedom. The Bayesian inverse problem is thus to characterize an implicitly defined distribution in a high-dimensional space. This is a computationally demanding problem that requires scalable and efficient numerical methods be used throughout: in discretizing the forward model; in solving the resulting nonlinear equations; in solving the Bayesian inverse problem; and in propagating the uncertainty encoded in the posterior distribution of the inverse problem forward onto important quantities of interest. To address discretization, a hybrid parallel adaptive mesh refinement format is designed and implemented for ice sheets that is suited to the large width-to-height aspect ratios of the polar ice sheets. An efficient solver for the nonlinear Stokes equations is designed for high-order, stable, mixed finite-element discretizations on these adaptively refined meshes. A Gaussian approximation of the posterior distribution of parameters is defined, whose mean and covariance can be efficiently and scalably computed using adjoint-based methods from PDE-constrained optimization. Using a low-rank approximation of the covariance of this distribution, the covariance of the parameter is pushed forward onto quantities of interest.