ABSTRACT. The so-called magnetoquasistatic approximation or eddy current problem is a standard approach to simulate the behavior of electrical machines. Thereby Maxwell's equations are simplified by neglecting the displacement current density in Ampère's circuit law. The resulting time-dependent eddy current problem is typically solved by a time-stepping approach, which sequentially solves one time step after the other. For long time periods, as for example at the startup of an electrical machine, the time-stepping approach leads to high computational costs. In contrast, parallel-in-time methods enable parallelism by computing multiple time steps simultaneously and, therefore, are promising techniques for reducing the simulation time of this long range time-dependent problems.

We apply two parallel-in-time algorithms (Parareal and Multigrid Reduction in Time (MGRIT)) to the eddy current problem of a two-dimensional coaxial cable model problem and compare the runtime results of the time-parallel solutions to the runtime of the time-stepping approach.

ABSTRACT. Transient PDE constrained optimization is challenging due to repeated forward and backward simulation. Using current approaches, when a simulation uses only a parallel spatial decomposition near the strong scaling limit, the time to solution for the optimization problem cannot be decreased by adding more computational resources. Addressing this issue requires algorithmic advancement in optimization algorithms and solution methods. This talk proposes a parallel-in-time multigrid preconditioner for solving transient KKT systems arising in an inexact SQP algorithm. Our composite-step SQP algorithm utilizes inexact iterative solution of "benign" KKT systems, corresponding to a sequence of strictly convex quadratic programs. Its inexactness-handling mechanisms ensure global and fast local convergence. The preconditioner used to solve the KKT system is critical to the efficiency of this algorithm. Our preconditioner explicitly exposes continuity in time constraints using a time domain decomposition technique. These constraints are relaxed and the coupled forward-adjoint system is solved using a Krylov method preconditioned with a multigrid in time algorithm. This talk will present the preconditioner and show results indicating convergence the scalability of the approach. Results demonstrating speedups assuming spatial parallelism is saturated for two dimensional nonlinear PDEs will be presented.

ABSTRACT. Fast tensor-product solvers have been proposed by Maday and Ronquist in 2007. The parallelization in time relies on a block-diagonalization. In that process, the choice of the time steps is crucial, as it deteriorates both the precision and the condition number of the linear system. In recent works with Martin Gander we proposed an optimization strategy for the choice of these parameters, applied to the heat equation, discretized in time with backward Euler. We propose here a general analysis for the theta-method. We rely on the computation of the eigenvectors for the full discretization matrix, and very precise truncation error estimates. This is a joint work with Martin Gander (Genève) et Juliette Ryan (ONERA).

ABSTRACT. We present recent progress towards constructing a parareal algorithm for solving the rotating shallow water equations. These equations are the simplest to exhibit behaviour relevant for atmospheric modelling and as such are often used as a test bed for new numerical algorithms. The parareal algorithm has two ingredients: a cheap integrator (the coarse propagator) and a more accurate integrator (the fine propagator). The convergence rate, and hence efficiency, of this scheme is strongly affected by the accuracy of the coarse propagator. Our approach follows that of [1] where it was shown that averaging the nonlinear terms over some fast oscillations includes the effects of near-resonances, essential for accuracy and hence convergence. A key component of this scheme is the computation of the exponential of the linear operator corresponding to the fast linear waves in the system. Here we use the rational approximation approach of [2] (REXI). This requires the solution of an elliptic problem for each term. The REXI method is highly parallelisable as each term can be computed separately, but the solution of each term must be efficient, scalable and fast. In this talk, we will describe strategies for solving the elliptic problem and present parallel scaling results for REXI applied to a compatible finite element model for the linear shallow water equations on the sphere.

[1] Haut, T. and Wingate, B. An asymptotic parallel-in-time method for highly oscillatory PDEs. SIAM Journal on Scientific Computing, 36.2 (2014), A693-A713

[2] Haut, T. et al A high-order time-parallel scheme for solving wave propagation problems via the direct construction of an approximate time- evolution operator. IMA Journal of Numerical Analysis, 36(2) (2015), 688-716

ABSTRACT. Residual neural networks (ResNets) are a type of deep neural network that exhibits excellent performance for many learning tasks, e.g., image classification and recognition. Mathematically, ResNet architectures are equivalent to a forward Euler discretization of a nonlinear initial value problem, where the time-dependent control variables are the network weights and each network layer is associated with a time-step. Thus, ResNet training can be represented as an optimal control problem of the associated dynamical system. However, both ResNet training and similar time-dependent optimal control problems from engineering applications often suffer from prohibitively long run-times because of the many sequential sweeps forwards and backwards across layers (i.e., time-steps) to carry out the optimization. One remedy for this issue is to use parallel-in-time methods, which have shown notable improvements in scalability for some engineering applications. This work demonstrates the use of one such parallel-in-time technique (multigrid-reduction-in-time) for the efficient and effective training of ResNets. The proposed algorithms replace the classical (sequential) forward and backward propagation through the network layers by a parallel nonlinear multigrid iteration applied to the layer domain. This adds a new dimension of parallelism across layers that is attractive when training very deep networks. From this basic idea, we derive multiple layer-parallel methods, which allow for the simultaneous training of all network layers. The most efficient version employs a simultaneous optimization approach based on inexact gradient information. Using supervised classification examples, we demonstrate that our approach achieves a similar training performance to traditional methods, but enables layer-parallelism and thus provides speedup over layer-serial methods.

ABSTRACT. With the increased interest in electronic vehicles, there is a focus on increasing the performance of Li-ion batteries. Two ways to increase performance are: increase the battery’s ability to charge faster or increase its charge density. Either method can be explored with the aid of highly accurate simulations.

A Li-ion battery contains three domains: two electrodes (anode and cathode) and the electrolyte that connects them. The geometry of these domains tend to be extremely complex. Solving the full scale system requires a fast iterative solver with an appropriate preconditioner.

The mathematical model describes the Lithium concentration and potential in all three domains at the microscopic scale. The transport of Lithium ions at the interface of an electrode and the electrolyte is governed by the nonlinear Butler-Volmer equation. In the electrodes, this system is very Poisson-like, and is a perfect candidate for an Algebraic Multigrid preconditioner (AMG). However, within the electrolyte, the Poisson-like system is augmented with a second-order, nonlinear, off-diagonal term that has the potential to ruin the AMG performance. The current best solution strategy is splitting the system into two blocks, one for the concentration and one for the potential. Then, perform block GMRES preconditioned with AMG. This solve is facilitated using the FEniCS and PETSc libraries.

The goal of this presentation is to discuss preconditioning methods for modeling Li-ion batteries. It will start by introducing the system of equations that describe the transport of Lithium Ions throughout the battery. Then, the performance of multiple solution strategies will be presented. Finally, we hope to ideas for other preconditioning strategies.

ABSTRACT. Neutron transport equations have been widely used in the study of nuclear reactor behaviors. With the advancements of supercomputers, parallel simulations of the neutron transport equations become popular, but it is a computationally challenging task since the scalable solver is difficult to design and develop for the neutron transport equations that is defined on a seven-dimensional space. In this work, we propose a scalable eigenvalue solver together with a subspace coarsening AMG method for the numerical solution of neutron transport equations. A nonlinear eigenvalue problem in the neutron transport equations is first reformed as a regular nonlinear system of equations, and then the nonlinear system is calculated using an inexact Jacobian-free Newton method. In order to speedup the simulation, an efficient preconditioner is applied during each Newton iteration for calculating the Jacobian system of equations. In the preconditioner, a subspace coarsening algorithm is introduced to reduce the preconditioner setup time, and the resulting subspace interpolation/restriction is extended to construct a full space AMG method for solving the coupled Jacobian system. We numerically show that the proposed algorithm is scalable with more than 10,000 processor cores for realistic reactor simulations with more than billions of unknowns.

ABSTRACT. Simulation of wave interactions with configurations containing extremely large numbers of individual particles are important in diverse applications including oceanography and atmospheric science applications. Simulations need to account for reflections between all particles, the number of which grows with the square of the number of particles. Consequently such simulations have been limited to only a few thousand particles. We present an efficient solver based on a Krylov subspace method combined with a fast algorithm for computing matrix vector products whose complexity (memory and CPU time) grows only linearly with the number of particles. Our matrix vector product scheme is based on a novel numerically stable expansion of the field radiated by each particle in spherical wavefunctions. We demonstrate our algorithm by simulating multiple scattering for configurations with more than a quarter of a million individual particles and associated dense complex linear systems with more than one hundred million unknowns.

ABSTRACT. We develop a fully-coupled, fully-implicit approach to phase field modeling of solidification for additively manufactured materials. Predictive simulation of solidification in pure metals and alloys remains a significant challenge in the field of materials science, as micro-structure formation during the solidification of a material plays an important role in the properties of the solid material. Our approach consists of a finite element spatial discretization of the fully-coupled nonlinear system of partial differential equations at the microscale, which is treated implicitly in time with a preconditioned Jacobian-free Newton-Krylov (JFNK) method. The approach allows timesteps larger than those restricted by the traditional explicit CFL limit and is algorithmically scalable and efficient due to an effective preconditioning strategy based on algebraic multigrid and block factorization. Preconditioner performance is analyzed in terms of parallel scalability and efficiency on heterogeneous architectures.

ABSTRACT. This talk introduces a discretization-accurate stopping criterion of symmetric iterative methods for solving systems of algebraic equations resulting from finite element approximation. The stopping criterion consists of the discretization and the algebraic error estimators, that are based on the respective duality error estimator and difference of two consecutive iterates. Iteration is terminated when the algebraic estimator is in the same magnitude as the discretization estimator. Numerical results for multigrid V(1,1)-cycle and symmetric Gauss-Seidel iterative methods are presented for the linear finite element approximation to the Poisson equations. A large computation reduction is observed comparing to the standard residual stopping criterion.

ABSTRACT. Hierarchical hybrid grids (HHG) have been proposed in conjunction with multigrid solvers to promote the efficient utilization of high performance modern computer architectures. While HHG meshes provide some flexibility for unstructured applications, most of the multigrid calculations can be accomplished using structured grid ideas and kernels. In this paper, we generalize the HHG idea in two ways to facilitate its adoption within the unstructured finite element community. First, we introduce a mathematical framework that can address HHG meshes as well as more general types of block structured meshes. The framework illustrates the relationship between the new partially structured MG algorithms and more traditional MG algorithms, and helps provide a guide for implementing new semi-structured multigrid algorithms. Overall, our implementation focuses on significantly reducing the overall effort required by application scientists to adapt their mature finite element capabilities to an HHG solver. This is accomplished by applying a dis-assembly process to a standard finite element matrix. Second, we modify the HHG solver so that it can address significantly more complex meshes. In particular, conventional HHG meshes are constructed by uniformly refining an unstructured mesh. Instead, we consider extensions such that the underlying finest mesh can contain a few completely unstructured regions. These extensions allow finite element practitioners to use unstructured meshes in specific areas (e.g., near jagged interfaces) where they are most natural.

ABSTRACT. Due to their grid-independent convergence multigrid algorithms are popular for the solution of elliptic PDE which have been discretised with (higher-order) Discontinuous Galerkin (DG) methods.

To fully utilise computational resources, efficient and algorithmically optimal implementations are necessary. Traditionally, a given problem is solved by assembling a system of sparse equations and inverting the resulting matrix algebraically. On modern manycore chip architectures with a poor FLOP-to-bandwidth ratio this approach is very expensive with computational complexity of O(N^2) where N denotes the number of unknowns per cell. Matrix-free sum factorisation techniques in d dimensions reduce the computational complexity from O(N^2)=O(n^{2d}) to O(d*n^{d+1}), where n is the number of unknowns in one direction; this leads to significant speedups in the operator application [Muething, Bastian,Piatkowski (2017), arXiv:1711.10885].

In multigrid smoothers it is also necessary to invert block-diagonal matrices. Matrices D_e of size n^d*n^d are assembled in each cell e and inverted exactly. This requires O(n^{2d}) bandwidth-bound operations and quickly becomes the bottleneck of the solver as the order n increases, obliterating gains from the matrix-free operator application. To circumvent this, we solve the system D_e*x=y iteratively by implementing multiplication with D_e in a matrix-free way. This reduces the cost from O(n^{2d}) to O(n_{iter}*d*n^{d+1} where n_{iter} is the number of iterations to solve the system in each cell.

We solve linear convection-diffusion systems and demonstrate the algorithmic and computational efficiency of the method for a hybrid multigrid algorithm with hp-coarsening, similar to [Bastian et al. (2012), Num. Lin. Alg. with Appl. 19 (2), pp. 367-388]: on the finest level a matrix-free block smoother is applied to the high-order system, and the low-order system on the coarser levels is solved with AMG. We demonstrate the efficiency of our EXADUNE based implementation for a range of elliptic PDEs, including the SPE10 benchmark for subsurface flow.

ABSTRACT. For large scale computing, memory size and the cost of accessing large memory repeatedly can easily become the dominating computational bottleneck. This is often true for efficient multigrid methods, which can ideally solve a partial differential equation in only a moderate number of work units. A work unit is here understood in the sense of Achi Brandt's textbook efficient paradigm as one application of the discretized operator. Efficient solvers must therefore attempt to solve a given problem in as few work units as possible, but they must also try to make a work unit itself as cheap as possible. Often the cost of a work unit can be reduced by matrix-free techniques. Here hierachical hybrid grid (HHG) structures can be used to design algorithms that permit a cheap re-calculation of the nonzero matrix elements whenever they are needed. The key to these techniques is to exploit locally structured sub grids in a globally unstructured mesh. We will discuss two techniques. One method uses a scaling of the matrix rows (aka stencils), another promising approach uses moderate order polynomials to approximate the matrix entries. In both cases the correct stiffness matrix is only approximated so that the corresponding error can be analyzed as a so-called variational crime in a finite-element setting. Both techniques are useful in massively parallel computing. Their use for large scale simulations will be demonstrated for applications arising in geophysics.

ABSTRACT. Getting good performance and good resource utilization on next generation platforms (NGPs) is a current challenge for linear solvers. In order to overcome this challenge solver code bases need to be re-factored to utilize more effective kernels. This opens the door for two approaches: first one can work on new implementation of current algorithms taking care of exposing more on node parallelism or new algorithms can be explored. In this talk I would like to demonstrate how partially structured grids that exist in many applications can be leveraged to tackle this challenge. New algorithms for structured grids have been implemented in MueLu the multigrid framework within the Trilinos project. Their performance (numerical and algorithmic) are compared to classic AMG methods typically used in MueLu.

ABSTRACT. Line and plane relaxation are important components for robust variational multigrid methods on structured grids when using standard coarsening with anisotropic problems. These methods constitute the majority of time spent in a multilevel solve--driving the performance of the solver. Tridiagonal systems used in line relaxation are solved using a memory efficient, multilevel parallel partitioning algorithm. In plane relaxation, 2D multilevel cycles are run on each plane. This results in a series of independent parallel cycles. To improve parallel efficiency, the uniform communication pattern of these cycles can be exploited to aggregate communication across planes on a process. A straightforward application of this approach, when used with coarse-grid redistribution at scale, is limited by the maximum number of communicators in MPI so coordination is needed among planes. In this talk, we will present an approach for automating this aggregation using lightweight user-level threads.

ABSTRACT. Magnetohydrodynamics (MHD) models describe a wide range of plasma physics applications, from thermonuclear fusion in tokamak reactors to astrophysical models. These models are characterized by a nonlinear system of partial differential equations in which the flow of the fluid strongly couples to the evolution of electromagnetic fields. As a result, the discrete linearized systems that arise in the numerical solution of these equations are generally difficult to solve, and require effective preconditioners to be developed. This talk investigates monolithic multigrid preconditioners for a one-fluid, viscoresistive MHD model in two dimensions that utilizes a second Lagrange multiplier added to Faraday's law to enforce the divergence-free constraint on the magnetic field. We consider the extension of a well-known relaxation scheme from the fluid dynamics literature, Vanka relaxation, to this formulation. To isolate the relaxation scheme from the rest of the multigrid method, we utilize inf-sup stable elements for both constraints in the system, and a geometric multigrid approach with finite-element interpolation operators and rediscretization for the coarse grid construction. An additive variant of the Vanka-scheme is discussed, and parallel numerical results are shown for the Hartmann flow problem, a standard steady-state test problem in MHD, as well as a time-dependent Island Coalescence problem.

ABSTRACT. The numerical simulation of flows associated with laser-induced phase change in metal additive manufacturing, supercritical fluids, and the slow cook-off of energetic materials present new challenges. Specifically, these flows require a fully-implicit, all-speed compressible formulation, since rapid density variations occur due to phase change over long hydrodynamic time-scales.

We investigate the preconditioning for a high-order all-speed compressible Navier-Stokes solver in the ALE3D multi-physics code that addresses such challenges. The equations are discretized with a reconstructed Discontinuous Galerkin (rDG) method and integrated in time with fully implicit discretization schemes. The resulting set of non-linear and linear equations are solved with a robust Jacobian-Free Newton-Krylov (JFNK) framework.

In the limit of large acoustic CFL number, there is a tight coupling between the low and high order degrees of freedom, leading to a 3x3 block system per equation for a 2nd-order scheme and a 6x6 block system per equation for a 3rd-order scheme. To robustly solve the large ill-conditioned block systems, we utilize a nested Schur complement approach to reduce the NxN high-order blocks to N 1x1 blocks. The nested Schur complement preconditioner is a polynomial multigrid (p-multigrid) reduction technique, where the reduced systems can be effectively solved with an algebraic multigrid method. Numerical results are shown for the p-multigrid Schur complement preconditioned Newton-FGMRES solver in the limit of large CFL numbers for compressible internally heated convection and a supercritical flow problem.

ABSTRACT. A new method for high-order discretization of generally discontinuous solutions in multifluid flow simulations will be presented. The method is based on the arbitrary-order Discontinuous Galerkin (DG) space discretization combined with the implicit method-of-lines time discretization. We use the DG-based level set to evolve multimaterial contacts. The level set is used to accurately integrate conservation laws, accounting for presence of multiple solution fields in each "mix" computational elements, using innovative curvilinear cutcell algorithm, without any assumption of mechanical or thermal equilibrium. Combining high-order least-squares based reconstruction with interfacial jump conditions, the reconstructed solution fields allows to implement very complex interfacial physics, accurately accounting for jumps in material properties (density, viscosity, thermal conductivity) and surface tension without the need to add numerical diffusion to regularize the solutions. The method is demonstrated to work robustly, with the DG space discretization up to the 4th-order. The method is based on the fully-implicit Newton-Krylov solution procedure, which was essential for robustness and enabling conservation of the mixture. Numerical examples cover a wide range of applications related to laser directed energy and fluid-structure interactions.

ABSTRACT. Simulation of multiphase poromechanics involves solving a multi-physics problem in which multiphase flow is tightly coupled with poromechanical deformation. To capture this dynamic interplay, fully implicit methods, also known as monolithic approaches, are usually preferred. The main bottleneck of this strategy is the cost of solving the linear systems resulting from discretization of the problem. Because of the strong coupling present in the continuous problem, efficient techniques such as algebraic multigrid (AMG) cannot be directly applied to the discrete linear systems. In this work, we present our efforts in developing an algebraic framework based on multigrid reduction that is suited for tightly coupled systems of PDEs. Using this framework, the decoupling between the equations is done algebraically through defining appropriate interpolation and restriction operators. One can then employ existing solvers for each of the decoupled blocks or design a new solver based on knowledge of the physics. We will demonstrate the applicability of our framework to multiphase flow coupled with poromechanics. We show that the framework is flexible to accommodate a wide range of scenarios, as well as efficient and scalable for large problems.

ABSTRACT. For problems in Numerical Weather Prediction (NWP), time to solution is a critical factor. Semi-implicit time-stepping methods can speed up geophysical fluid dynamics simulations by taking larger time-steps than explicit methods. This is possible because they treat the fast (but physically less important) waves implicitly, and the time-step size is not restricted by the CFL condition for these waves. One disadvantage of this method is that an expensive linear solve must be performed at every time step, however, using an effective preconditioner for an iterative method significantly reduces the computational cost of this solve, making a semi-implicit scheme faster overall.

Higher-order Discontinuous Galerkin (DG) methods are known for having high arithmetic intensity making them well suited for modern HPC hardware, but are difficult to precondition due to the large number of coupled degrees of freedom. This coupling arises since the numerical flux introduces off diagonal artificial diffusion terms. By using a hybridised DG method we can eliminate the original coupling and instead couple the equations to a smaller global system on the trace space, which is easier to precondition. This is achieved by considering the numerical flux variables which only lie on the facets of the mesh. Recent work by Kang, Giraldo and Bui-Thanh[1] solves the resultant system directly. However, this becomes impractical for high resolution simulations. We build on this work by solving the resultant system using a non-nested geometric multigrid technique instead[2].

We discretise and solve the non-linear shallow water equations, an important model system in geophysical fluid dynamics, and demonstrate the effectiveness of the multigrid preconditioner for a semi-implicit IMEX time-stepper. The method is implemented in the SLATE language, which is part of the Firedrake project. Firedrake is a Python framework for solving finite element problems via code generation.

References:

[1] Kang, Shinhoo and Giraldo, Francis X and Bui-Thanh, Tan. IMEX HDG-DG: a coupled implicit hybridized discontinuous Galerkin (HDG) and explicit discontinuous Galerkin (DG) approach for shallow water systems arXiv preprint arXiv:1711.02751, 2017

[2] Cockburn, Bernardo and Dubois, Olivier and Gopalakrishnan, Jay and Tan, Shuguang. Multigrid for an HDG method IMA Journal of Numerical Analysis 34(4):1386–1425, 2014

ABSTRACT. An efficient hp-multigrid scheme is presented for local discontinuous Galerkin (LDG) discretizations of elliptic problems, formulated around the idea of separately coarsening the underlying discrete gradient and divergence operators. We show that traditional multigrid coarsening of the primal formulation leads to poor and suboptimal multigrid performance, whereas coarsening of the flux formulation leads to optimal convergence and is equivalent to a purely geometric multigrid method. The resulting operator-coarsening schemes do not require the entire mesh hierarchy to be explicitly built, thereby obviating the need to compute quadrature rules, lifting operators, and other mesh-related quantities on coarse meshes. We show that good multigrid convergence rates are achieved in a variety of numerical tests on uniform Cartesian grids as well as for curved domains using implicitly defined meshes and for multi-phase elliptic interface problems with complex geometry.

ABSTRACT. In this work we discuss robust linear solvers of a stabilized discretization of the poroelastic equations. The discretization is well-posed with respect to the physical and discretization parameters, and thus provides a framework to develop block preconditioners that are robust with respect to such parameters as well. We construct these preconditioners for the stabilized discretization and a perturbation of the stabilized discretization that leads to a smaller overall problem. Numerical results confirm the robustness of the block preconditioners. Time permitting, we also discuss a monolithic geometric multigrid method for solving the stablized discretization and compare its performance with the block preconditioners.

ABSTRACT. In many models used in engineering and science, material properties are uncertain or spatially varying. For example, in geophysics, and porous media in particular, the uncertain permeability of the material is modelled as a random field. Depending on the material, these random fields can be highly anisotropic. Efficient solvers, such as the Multiple Semi-Coarsened Multigrid (MSG) method, see [11, 12, 13], are required to compute solutions for various realisations of the uncertain material. The MSG method is an extension of the classic Multigrid method that uses additional coarse grids that are coarsened in only a single coordinate direction. In this sense, it closely resembles the extension of Multilevel Monte Carlo (MLMC), see [4], to Multi-Index Monte Carlo (MIMC), see [7]. We present an unbiased MIMC method that reuses the MSG coarse samples, similar to the work in [9] and [16]. Our formulation of the estimator can be interpreted as the problem of learning the unknown distribution of the number of samples across all indices, and, in this sense, unifies the previous work on adaptive MIMC from [15] and unbiased estimation from [14]. We analyse the cost of this new estimator and present numerical experiments with various anisotropic random fields, where the unknown coefficients in the covariance model are considered as hyperparameters. We illustrate its robustness and superiority over standard MIMC without sample reuse.