ABSTRACT. Scalable solution of time-harmonic high-frequency wave propagation problems—including acoustic Helmholtz, elastic Helmholtz, and time-harmonic Maxwell—remains a challenge in mathematics and scientific computing. One difficulty is that classical discretization techniques (e.g., Galerkin finite elements, finite difference, etc.) yield indefinite discrete systems that preclude use of many scalable solution algorithms. Significant progress has been made to develop specialized preconditioners for high-frequency wave propagation problems, but robust and scalable solvers for general problems including non-homogenous media and complex geometries remain elusive. An alternative approach is to use minimum residual discretization methods—that yield definite discrete systems—and may be amenable to more standard preconditioners. The discontinuous Petrov-Galerkin (DPG) finite element methodology is one such minimum residual method with additional properties including mesh-independent stability and a built-in error indicator.

The present work details a scalable multigrid solver for general DPG systems (discretized with the standard energy spaces)—including high-frequency wave propagation and multiphysics problems. The DPG multigrid (DPG-MG) solver leverages the mesh-independent stability and built-in error indicator of the DPG methodology to define a hierarchy of hp-adaptive meshes on which multilevel preconditioning is performed. Mixed meshes with elements of all types (hexahedra, tetrahedra, prisms, and pyramids) as well as isotropic and anisotropic refinements are supported. A hybrid OpenMP/MPI implementation reveals that the DPG-MG solver is highly competitive for high-frequency wave propagation problems in both pre-asymptotic and asymptotic regimes. In particular, we demonstrate hp-adaptive solution of various three-dimensional time-harmonic wave propagation problems, including problems with more than 100 wavelengths, multiple billion degrees of freedom, non-homogeneous media, and general geometries (e.g. tokamak).

ABSTRACT. Algebraic multigrid (AMG) is a popular and effective solver for sparse linear systems arising from discretized PDEs. While traditional AMG methods are effective for a wide range of SPD systems, such methods often struggle for the nonsymmetric case. A recent AMG method, Approximate Ideal Restriction (AIR), addresses this issue by appealing to multigrid reduction ideas in order to effectively solve many nonsymmetric systems, often arising from problems with strong advection and anisotropy. While AIR is an effective method for many of these systems, AIR still has some drawbacks which we seek to address. We explore (i) reducing the possibly large operator complexity through faster coarsening, (ii) improving performance for strongly elliptic problems through consideration of the near kernel, and (iii) using local approximate inverses when constructing the restriction operator for efficiency. The goal is a more general and robust AMG method for symmetric and nonsymmetric systems.

ABSTRACT. Algebraic multigrid (AMG) method is proven to be an effective solver for symmetric positive definite (SPD) linear systems resulting from the discretization of many elliptic and parabolic PDEs. However existing AMG methods and heuristics mostly rely on the discretization matrix $A$ being symmetric and for this reason, highly nonsymmetric linear systems, which often arise from discretizing hyperbolic PDEs, remain a challenging problem for classical AMG methods and its variants. AMG for such problems is typically hard because the traditional coarsening, interpolation and relaxation strategies designed for symmetric positive definite systems can break down. In this work, we investigate reduction-based AMG based on a local approximation to the ideal restriction operator (AIR), as well as the bootstrap framework based on a multigrid eigensolver targeting nonsymmetric problems. As an application problem, we consider a time dependent relativistic drift-kinetic Fokker-Planck model in the phase space implemented with adaptive mesh refinement. We demonstrate how applying reduction based AMG can in some cases lead to improved convergence over classical AMG. The ultimate goal is to develop an algebraic multigrid solver which addresses various nonsymmetric discretizations, including time-dependent and steady-state advection-diffusion, ranging from the purely advective to the purely diffusive.

ABSTRACT. High-resolution, continental scale hydrologic models have increased in use and application. These models analyze interactions between hydrologic phenomena beyond catchment boundaries and across scales. A common problem for users of these types of models is that they are large and complex, with immense code bases that are difficult to navigate, especially for domain scientists with limited computer science backgrounds. We present a hydrologist’s process of debugging, running, and performance testing the integrated hydrologic model ParFlow-CLM at the continental scale. ParFlow is a physics-based, integrated hydrological model that simulates variably saturated subsurface and surface flow equations simultaneously. In this use case, ParFlow is coupled to the Common Land Model (CLM). The code is parallelizable and has been run on various high performance computing systems, making modeling at large spatial scales feasible. The continental scale application of ParFlow-CLM, ParFlow-CONUS version 2 (PF-CONUSv2), represents a domain that incorporates the entirety of the continental US and a 10 layer deep subsurface, resulting in 78.5 million total unknowns. PF-CONUSv2 is being run on the NCAR Cheyenne supercomputer but we are experiencing slow simulation speeds that make research simulations inefficient to run. To measure and attempt to optimize performance we have conducted parallel scaling tests to determine the most efficient processor topology, as well as utilized a GNU profiling tool to better understand which components of the code may be slowing the program execution. While these processes might be straightforward to a computer scientist, the particulars of large multi-gridded models present their own unique challenges to be overcome. We hope this example application of multigrid methods can illuminate some workflows and problems domain scientists may experience.

ABSTRACT. Deep neural networks (DNNs), a special form of machine learning, are a powerful tool for many learning tasks, such as image recognition. Although significant progress has been made to speed up the training process of these DNNs, the training process often remains a bottleneck due to the increasing complexity and depth of the models. Mathematically, the training process can be viewed as an evolutionary equation where the weights of a network layer are adjusted at each step. A relatively new approach to accelerate the training process through parallelization is to apply the time-parallel MGRIT method to this equation, allowing simultaneous training of all network layers.

In this talk, we present the Torchbraid framework, which is built on PyTorch and interfaces the MGRIT library XBRAID for layer-parallel neural network training. We show various runtime results of layer parallelism on CPUs and GPUs, results of combining layer and data parallelism, and other features of the framework.

ABSTRACT. This talk will investigate optimized parameters for domain decomposition methods (DDMs) with multiple levels. Prior work has considered optimizing the parameters for domain decomposition methods (DDMs) in the one-level setting or in special cases such as structured-grid discretizations with regular subdomain construction. We extend this work to multilevel DDM by proposing multigrid graph neural networks (MG-GNN), a novel GNN architecture for learning optimized parameters in multilevel Optimized Restricted Additive Schwarz (ORAS). We train MG-GNN using a new unsupervised loss function, enabling effective training on small problems that yields robust performance on unstructured grids that are orders of magnitude larger than those in the training set. We show that MG-GNN outperforms popular hierarchical graph network architectures for this optimization and that our proposed loss function is necessary to achieve this improved performance.

ABSTRACT. Solving large scale linear systems arising from the discretization of partial differential equations is a difficult and computationally expensive task, due to the large amount of degrees of freedom and the large condition number. The recently developed Hybrid Iterative Numerical Transferrable Solver (HINTS) provides an efficient solution strategy by combining a standard iterative solver (smoothing high frequency errors) with a Deep Operator Network (reducing low frequency errors). In this work we discuss various strategies for using the HINTS as a preconditioner for enhancing convergence of Krylov methods. We demonstrate the efficiency and robustness of the proposed preconditioner for several numerical examples, such as the non-SPD Helmholtz problem, Darcy flow, etc. We provide analysis supporting the results and show that the proposed preconditioners compete even with the state-of-the-art algebraic Multigrid preconditioners.

ABSTRACT. Patch-based multigrid smoothers are a class of smoothers for multigrid methods which solve local dense problems and reconstruct an approximate solution for the global problem. These smoothers are frequently used for linear systems obtained from high-order discretizations and multiphysics applications, as they contain dense local couplings which are not effectively handled by traditional smoothers. While patch-based smoothers are generally more accurate than traditional methods for such problems, they are often more expensive. We present data-driven approaches for compressing and constructing reduced representations of patches. These reduced representations greatly increase the smoother speed with little or no increase in iteration count. We present a wide variety of numerical experiments on various applications to verify our method.

ABSTRACT. We present a class of multi-level algorithms for unconstrained nonlinear optimization, which does not require the evaluation of the objective function [1]. The choice of avoiding the evaluation of the objective function is intended to make the algorithms less sensitive to noise, while the multi-level feature aims at reducing their computational cost. The proposed algorithmic class encompasses several novel methods as well as the well-known AdaGrad method as a specific (single-level) instance. The global convergence rate of the proposed algorithms is analyzed and their numerical behavior in the presence of noise is illustrated in the context of training deep neural networks for supervised learning applications. Comparison with the state-of-the-art optimizers is also performed, demonstrating encouraging results.

[1] S. Gratton, A. Kopanicakova, Ph. L. Toint. Multilevel objective-function-free optimization with an application to neural networks training. Under review (2023)

ABSTRACT. We consider the application of the multigrid reduction-in-time (MGRIT) algorithm to the acoustic equations, a first-order, linear hyperbolic system. It is well known that solving hyperbolic problems with MGRIT is challenging, with the key to success being to devise a coarse-grid problem that approximates closely enough the fine-grid problem. In this work, we develop such a coarse-grid problem by using a semi-Lagrangian-like discretization which applies to the characteristic variables of the system. This coarse-grid discretization involves a correction term that ensures its truncation error mimics closely enough that of the fine-grid discretization. This extends recent work of ours in which we applied related ideas to solve scalar linear advection equations with MGRIT.

ABSTRACT. We describe a new strategy for solving highly oscillatory two-dimensional Vlasov-Poisson systems by means of a specific version of the parareal algorithm. We use reduced models, obtained from the two-scale convergence theory, for the coarse solving. These models are useful to approximate the original Vlasov-Poisson model at a low computational cost since they are free of high oscillations. The models are numerically solved with a particle-in-cell method. We provide a thorough analysis of the efficiency of the parareal algorithm on the basis of the ideal speedup. Next, within a distributed memory paradigm we discuss the performance of the parallel processing over the time slices.

ABSTRACT. Kinetic plasma simulations face the challenge of dealing with multiple time-scales that vary widely, resulting in stability limitations and the need to take numerous time-steps. Although parallel-in-time (PIT) methods are attractive, the standard PIT approach involves spatiotemporal coarsening, which can lead to iterations that are far removed from the physical solution or even unstable. This hinders convergence and limits speed-ups. To overcome this issue, we propose a novel PIT strategy that employs a coarsening of models, rather than meshes. This approach accelerates the time-to-solution for high-dimensional kinetic plasma models by using their much cheaper, lower-dimensional fluid approximations. We have tested this method on the two-stream instability and obtained promising initial results.

ABSTRACT. This work presents the first steps towards fast topology optimisation of time-dependent problems with many timesteps.

Topology optimisation is a process for automatically generating design concepts for complex engineering problems using a combination of simulations and optimisation algorithms. It is an iterative process, where the performance of the system is analysed using one or several simulations. The design is then updated until the performance is satisfactory, a demanding process usually requiring around 100-1000 simulations. This means that each simulation must be as fast as possible in order to reduce the overall time-to-solution for the user.

The current state-of-the-art is somewhat limited to steady-state (static) problems, since time-dependent (transient) simulations generally take a very long time due to the inherent serial nature of stepping through time to resolve the system behaviour. Practically, the computational time of time-stepping methods is infeasible in both industrial and academic timescales, because even simple time-dependent problems can take days to optimise. To overcome this limitation, this work considers a space-time formulation and solves the full time-dependent system using an "all-at-once" approach. Parallelisation and multigrid in space-time is applied to reduce the time-to-solution using distributed-memory computing clusters.

In this work, a continuous Galerkin finite element formulation with piecewise-linear basis functions is used to discretise the thermal diffusion equation with varying coefficients. Simple space-time "geometric" coarsening is used in conjunction with a Galerkin-projection geometric multigrid preconditioner. It is shown that a simple, and somewhat naive, implementation of the space-time framework for two-dimensional (spatial) problems leads to a 10% time-to-solution compared to time-stepping, albeit at a cost of 6 times the computational resources (core-hours). The simple setup works very well despite the issue of large contrast in coefficients in space (topology optimisation), even for non-linear problems with large contrast also in time (phase change).

ABSTRACT. Optimal control problems governed by time-dependent partial differential equations (PDEs) lead to large-scale discrete optimization problems. The size of these problems is proportional to the size of the discrete time domain. One approach to solve such problems is the use of algorithms based on the augmented systems. In the context of time-dependent problems, an augmented system is the KKT system for a convex linear-quadratic problem. However, in most numerical approaches the time discretization serializes the solution process and introduces a bottleneck. In the strong scaling limit, this bottleneck cannot be overcome by additional parallelization in space.

To accelerate the solution of linear-quadratic optimal control problems governed by PDEs, we propose a time-domain decomposition approach that introduces time parallelism into the optimization algorithm. This is achieved by introducing auxiliary state variables for each time interval and impose time-continuity constraint. Additionally, auxiliary variables are incorporated into the objective function. We demonstrate the potential effectiveness of our scheme as a time-parallel preconditioner for linear systems arising in inexact SQP.

ABSTRACT. In nuclear fusion reactors, neutral particles play an important role in shielding the reactor walls from the hot plasma in which the fusion reactions take place. The resulting coupled physics problem is high dimensional and has a multiscale nature. At the microscopic scale, the particle behaviour is captured in detail by kinetic equations, which can be solved using stochastic Monte Carlo methods. At the macroscopic scale, a limited number of macrosopic quantities of interest can be modeled using approximate fluid models, which can be solved using finite volume methods. In high collisional regimes, the microscopic scale description becomes increasingly more expensive to simulate. At the same time, the macroscopic scale description becomes increasingly more accurate. In this talk, we combine these two observations in a new two-phases hybrid method for the treatment of the neutral particles.

ABSTRACT. In many applications of interest, strongly anisotropic diffusion operators have to be solved with a high degree of accuracy in order to avoid excessive pollution across the direction of anisotropy. One way of doing so is by formulating the operator as a mixed system, with an auxilary variable denoting the diffused field's derivative in the direction of anisotropy.

The highly anisotropic nature of the problem makes it difficult to solve with standard multigrid methods. Here we present a novel solver technique for resolving this issue for problems with open field lines. By exploiting the fact that the auxiliary variable setup splits the diffusion formulation into two transport operators, we apply AIR - an algebraic multigrid solver well adapted for transport problems - to the mixed formulation. Here we highlight the formulation of the anisotropic diffusion operator, identify the solver technique and demonstrate the solver's effectiveness on some test problems implemented in Firedrake.

ABSTRACT. The robust, scalable simulation of flowing electrochemical systems is increasingly important due to the synergy between intermittent renewable energy and electrochemical devices. The high Peclet regime of many such applications prevents the use of off-the-shelf discretization methods. In this presentation, we propose a high-order DG scheme for the electroneutral Nernst-Planck equations. The chosen charge conservation formulation allows for appropriate penalization in the DG scheme and enables different treatments in the preconditioner: AMG for the potential block and ILU-based methods for the advection-dominated concentration blocks. For higher-order cases, a p-multigrid method is used for the potential block. Strong scaling results for different preconditioning approaches are shown for a large 3D flow-plate reactor example.

Work performed for U.S. DOE by LLNL under contract DE-AC52-07NA27344, supported by LLNL-LDRD projects 19-ERD-035, 22-SI-006. LLNL-ABS-844505

ABSTRACT. In this talk I will discuss modeling fluid flow through a deformable porous medium. I will discuss an existing approach based on Biot's consolidation model. Then I will introduce a system derived using energetic variational approach and discuss numerical methods and simulations.

ABSTRACT. We introduce a new multigrid (MG) method for solving variational inequality problems. Our proposed MG method extends the monotone MG method to also handle the generic linear bound constraints. First, we will discuss a method to decouple linear constraints by projecting them onto a new basis. To handle these decoupled constraints, we introduce a modified projected Gauss-Seidel method as smoother in our MG method. This method allows us to enforce the linear constraints, identify the active set, and monotonically minimize the energy function. We will discuss the necessary modification to the standard transfer operators. The global convergence of our MG method is ensured by combining the modified smoother and the transfer operators. We test our MG method for contact problems arising from the unfitted finite element discretization. Through a series of numerical examples, we showcase the robustness and efficiency of our proposed MG method and its level-independent convergence property.

ABSTRACT. Applications such as fault-tolerant algorithms favour multigrid methods for local recovery owing to their convergence properties. However, a family of shrinking (from coarse to fine) local meshes may present in some situations, such as performing parallel multigrid V-scheme on adaptively refined meshes. This local mesh structure results in convergence deterioration. We derive a simple adjustment on the standard inter-mesh transfer operators (natural injection and its adjoint) in one-dimension. Convergence improvement is observed in numerical experiments for the new set of operators.

ABSTRACT. This work is on efficient and fast reduced basis method for the computing the action of rational approximation of matrix functions. A direct computation of the action of such a rational approximation requires solving multiple systems with large scale sparse matrices. Such problems are typically of the form $(A + d_j M)u = f$, $j=1:n_p=\mbox{number(poles)}$, $0<d_1<d_2<\ldots<d_{n_p}$, where $A$ is a stiffness matrix and $M$ is a mass matrix.

Our method constructs the reduced basis from the first few directions obtained from the Preconditioned Conjugate Gradient applied to one of these linear systems. As we observe, only a small number of directions (5-20) are needed to approximate all of remaining inverses. This reduces the cost dramatically because: (1) We only use two of the large scale problems to construct the basis; and (2) the rest of the problems are restricted to this subspace and they are much smaller in size. We test our algorithms in problems of the form $A^s u=f$.

ABSTRACT. Local Fourier analysis (LFA) is an invaluable tool to assess the convergence behavior of smoothing iterations and multigrid methods on structured domains. In this talk we discuss the integer linear algebra foundation of LFA and more generally discrete Fourier transforms (DFT) on structured, discrete and finite domains in general form, i.e., lattice and crystal tori. We disseminate the definition of a canonical orthogonal basis of Fourier wave vectors and the intimate connection between position and frequency spaces.

As a result we show that LFA is merely an exercise of mapping the introduced structural foundations of DFT to shift invariant linear operators, i.e., multiplication operators/convolutions on crystalline structures. Exploiting this connection we give an algorithmic and thus fully automatic description of LFA relying on (not necessarily matching) definitions of the involved multiplication operators.

ABSTRACT. We study the combination of proximal gradient descent with multigrid for solving a class of possibly nonsmooth strongly convex optimization problems. We propose a multigrid proximal gradient method called MGProx, which accelerates the proximal gradient method by multigrid, based on utilizing hierarchical information of the optimization problem. MGProx applies a newly introduced adaptive restriction operator to simplify the Minkowski sum of subdifferentials of the nondifferentiable objective function across different levels. We provide a theoretical characterization of MGProx. First we show that variables at all levels exhibit a fixed-point property at convergence. Next, we show that the coarse correction is a descent direction for the fine variable in the general nonsmooth case. Lastly, under some mild assumptions we provide the convergence rate for the algorithm. In the numerical experiments, we show that MGProx has a significantly faster convergence speed than proximal gradient descent and proximal gradient descent with Nesterov's acceleration on nonsmooth convex optimization problems such as the Elastic Obstacle Problem.