ABSTRACT. Multigrid Reduction in Time (MGRIT) is one of the parallel-in-time methods and uses coarse-grid time integrators with the rediscretization approach. This approach works well for parabolic PDEs but leads to poor convergence of MGRIT for hyperbolic PDEs. The coarse-grid operator optimization proposed by De Sterck et al. shows that MGRIT works efficiently even for hyperbolic PDEs by constructing coarse-grid time integrators appropriately. Based on this idea, we propose a new construction method for coarse-grid time integrators using parameterized Runge-Kutta methods, which regard the coefficients of the Butcher tableau as parameters and use multistage zeroth-order schemes dedicated to MGRIT convergence. The convergence improvement is confirmed in three aspects: convergence rates, eigenvalue distributions, and stability regions. Numerical experiments show speedup results for time-stepping methods combined with spatial redistribution techniques that accelerate the coarse-level solvers.

ABSTRACT. In nonlinear data assimilation (DA) applications, iterative ensemble-based smoothers are advantageous because they incorporate the prediction model in the analysis. However, finding the optimal estimate can be expensive because it requires running the forward model multiple times over the DA window. We present a novel approach to reduce the cost by using multigrid-in-time (MGinT) within the iterative procedure of the maximum likelihood ensemble smoother (MLES). By performing the minimization process consecutively from coarse to fine grids over the DA window, we efficiently achieve an initial state for the finest-grid optimization. We demonstrate the implementation using the Lorenz-96 model and the Institute of Electrical and Electronics Engineers (IEEE) 39-bus test case. Our results show the potential of MGinT, while highlighting its capabilities and limitations. This will give insight to applying MGinT to DA methods to future problems with high dimensionality and strong nonlinearity.

ABSTRACT. For time-dependent problems, Parallel-in-Time (PinT) algorithms allow us to parallelize problems in the time dimension when space parallelization alone creates communication bottlenecks. Parareal and Multigrid Reduction-in-Time (MGRIT) are two examples of such PinT algorithms based on multigrid techniques, but they are not truly scalable since they coarsen the problem only in the time dimension.

We will focus on a more intrusive method: the Space-Time Multigrid algorithm with block-Jacobi relaxation introduced by Gander and Neumüller. This algorithm provides excellent scalability for parabolic problems up to millions of cores, while still being equally as fast as forward substitution on one core only.

We will show that the performance of this algorithm can be further improved by the optimization of the smoothing parameters. This will allow the algorithm to be up to twice as fast as the original one. Results will be presented for the heat equation discretized with Backward Euler.

ABSTRACT. We investigate optimized Schwarz domain decomposition methods in time for the control of transport equation. In the case of an internal control all over the domain, the optimization problem can be transformed into a system of two coupled PDEs. We then apply the time-domain decomposition strategy on this PDE system. Under Fourier analysis, we find optimized parameters for both one-sided and two-sided cases, in the latter case we get a much better convergence factor. We illustrate our results by numerical tests.

ABSTRACT. We present a parallel-in-time multigrid method for the solution of time-periodic parabolic evolution problems, with possibly nonlinear coefficients. After discretizing with finite elements in space and with backward Euler in time, we solve the resulting linear block-system of equations by means of a multigrid in time method. As a smoother we use a block-Jacobi method, which can be applied in parallel. We present numerical experiments based on applications in magneto-quasistatics.

ABSTRACT. The use of Machine Learning for scientific computations is accelerating rapidly. Innovative methods have been proposed to replace iterative solvers or enhance available solvers using learning. Based on the recent advances in operator regression, we propose a hybrid, iterative, numerical, and transferable solver (HINTS). HINTS combines standard relaxation methods and the Deep Operator Network (DeepONet). HINTS can provide faster solutions for a wide class of differential equations, while preserving the accuracy close to machine zero. HINTS can be applied to Multigrid to further enhance its efficiency and robustness. Through eigenmode analysis, we find that the individual solvers in HINTS target distinct regions in the spectrum of eigenmodes, resulting in a uniform convergence rate and hence exceptional performance of the HINTS overall. Moreover, HINTS applies to equations in multidimensions, is flexible with regards to computational domain, and transferable to different discretizations.

ABSTRACT. The discovery of fast numerical solvers is rapidly growing. New techniques have been introduced for many applications, especially in computational mechanics. Most numerical solvers are highly dependent on the problem geometry and discretization. The newly developed Hybrid Iterative Numerical Transferable Solver (HINTS) combines a standard solver with a neural operator to achieve better performance. In this work, we explore the "T" in HINTS, i.e., its geometry transferability properties. We first propose to directly employ HINTS built for a specific geometry to a different but related one without any adjustments. In addition, we propose using transfer learning to improve the convergence of HINTS. We conduct numerical experiments for Darcy flow and plane-strain elasticity problems. The results show that both the direct application of HINTS and the transfer learning enhanced HINTS can accurately solve these problems on different geometries, with a significant edge for the latter.

ABSTRACT. We consider the problem of how efficiently deep ReLU neural networks can approximate Sobolev functions. Existing sharp results are only available for the Sobolev spaces $W^s(L_\infty)$ in which the derivatives are measured in $L_\infty$. However, the solutions to PDEs can often only be bounded in Sobolev spaces $W^s(L_q)$ where the derivates are measured in a weaker $L_q$-norm with $q < \infty$. We study the optimal $L_p$-approximation rates for deep ReLU neural networks on the Sobolev class $W^s(L_q)$ for all $p,q$ and $s$. A necessary condition for approximation to be possible is that the Sobolev embedding condition holds. When the Sobolev embedding condition is strictly satisfied, i.e. when we have a compact embedding, our results give the optimal approximation rates for deep ReLU networks.

ABSTRACT. Successful deep neural networks for image classification consists of millions of learnable parameter, but are overparameterized. Consider the weight count as a function of the spatial dimensions of the input and the number of channels and layers. Using convolutions, the weight count scales linear in the spatial dimensions, but quadratic with the number of channels. Recent research employing multigrid concepts in deep neural networks has shown that weights can be saved significantly by appropriate weight sharing and a hierarchical structure in the channel dimension can improve the weight complexity to linear. Utilizing these findings, we introduce an architecture that establishes multigrid structures in all relevant dimensions, contributing a drastically improved accuracy-parameter trade-off. Experiments show that our structured reduction in weight count reduces overparameterization, in addition to improving accuracy over ResNet architectures on typical image classification benchmarks.

ABSTRACT. The use of Sparse Approximate Inverses (SPAI) has enjoyed success as a relaxation scheme for algebraic multigrid solvers, as when given proper parameters it can be robust and convergent when compared to other iterative methods. Traditional SPAI algorithms, however, optimize for a generic sparse inverse that do not take into account their usage; here, we investigate the use of gradient-based optimizers to obtain approximate inverses that specifically complement the use of multigrid solvers, as well as online optimization of the approximate inverse during the solve phase.

ABSTRACT. Solving optimization problems with transient PDE-constraints is computationally costly due to the number of nonlinear iterations and the cost of solving large-scale KKT matrices. These matrices scale with the size of the spatial discretization times the number of time steps. We propose a new $2$-level domain decomposition preconditioner to solve these linear systems when constrained by the heat equation. Our approach leverages the observation that the Schur-complement is elliptic in time, and thus amenable to classical domain decomposition methods. Further, the application of the preconditioner uses existing time integration routines to facilitate implementation and maximize software reuse. The performance of the preconditioner is examined in an empirical study demonstrating the approach is scalable with respect to the number of time steps and subdomains.

ABSTRACT. Optimal-control problems can give rise to boundary value problems (BVPs) whose solution contains the required control term. ParaDiag and ParaOpt are two algorithms to solve these optimality BVPs in a time-parallel way. In this talk, we generalize ParaDiag and apply its underlying preconditioner to improve ParaOpt's scalability. We take a broad view at the problem and formulate the methods applied to linear problems in the same framework. We show how, using this approach, convergence for both algorithms can be studied by finding eigenvalues of related matrices. By keeping the analysis generic, we were able to improve the state-of-the-art coarse propagator for ParaOpt and to directly apply ParaDiag's analysis to the ParaOpt preconditioners we developed. We illustrate the obtained analysis with numerical results.

ABSTRACT. In this talk, we present the convergence analysis of Paraopt algorithm applied to a linear optimal control problem. The Paraopt algorithm is a parallel in time method for solving a optimality system arising in partial differential equations (PDEs) constrained optimization. This solving required two time integration methods and depends in crucial way on a quadrature formula and a Runge-Kutta method. We use an operator norm analysis to show that the convergence rate depends on the discretization parameters and the smaller of the orders of the time integration methods. We illustrate this by some numerical examples.

ABSTRACT. In order to have better control of the dynamics of energy networks, an accurate simulation is mandatory. Since modern power systems are more and more complex, simulation methods have to adapt. Converters in a power system use high frequency switching to react in every changes, and simulating this requires small step size to resolve this behavior. There are already techniques to take switching with repeating, predictable patterns into account, i.e., approximating a pulse-width modulation (PWM) signal with a sine wave. In a power system there are also components whose switching behavior depends on the dynamics itself. In this presentation, the switch estimator (SE) applied to spectral deferred corrections (SDC) is presented together with examples to show its performance. After few iterations the SE predicts the time of the switch using interpolation and root finding, and time step is adapted. After restarting SDC with new step size SDC is able to resolve the singularity more accurate.

ABSTRACT. Space-time discretization methods are well suited to handle moving domains in electric machines such as the electric motor. The space-time domain is discretized at once, and in contrast to time stepping methods the movement can be captured by a fixed space-time mesh. For the solution of the magneto quasi-static Maxwell equations in the electric motor we formulate space-time finite element methods, considering a two-dimensional spatial domain and the time as the third dimension of the domain. As in the case of a fixed domain we are able to prove an inf-sup stability condition to ensure unique solvability. These space-time finite element methods allow for a direct as well as an iterative parallel computation for the magnetic flux density. Accordingly, the torque and the iron losses can be computed for the electric motor.

ABSTRACT. Numerical climate and weather prediction models require fast solvers for the linear equations that arise from semi-implicit timestepping. Many models use finite difference discretisations, for which the construction of multigrid preconditioners is straightforward. However, compatible finite element methods for the atmospheric equations of motion have recently attracted considerable interest. Preconditioning the saddle point system that arises in this approach is challenging since the velocity mass matrix is not diagonal. Hybridisable discretisations overcome this issue: weakly enforcing continuity of the velocity with Lagrange multipliers leads to a sparse system of equations, which has a similar structure as the pressure Schur complement. We describe how this hybridised system can be preconditioned with non-nested multigrid. Our approach requires a small number of iterations, shows excellent performance and scales to large core counts in the UK next-generation forecast model.

ABSTRACT. Irksome is a library based on the Unified Form Language (UFL) that enables automated generation of Runge--Kutta methods for time-stepping with finite element spatial discretizations of partial differential equations (PDE). Allowing users to express semidiscrete forms of PDE, it generates UFL for the stage-coupled variational problems to be solved at each time step. The Firedrake package then generates efficient code for evaluating these variational problems and allows users a wide range of options to deploy efficient algebraic solvers in PETSc. We will survey many of the features of Irksome, and then address the major challenge facing practitioners of fully implicit RK methods -- solution of the stage-coupled algebraic system. Many strategies that scale well with respect to mesh parameters and the number of stages are being explored in the literature, and we will present a "monolithic" framework to stage-coupled multigrid that appears quite promising. As an example, we will demonstrate this preconditioning strategy for RadauIIA-based time stepping of the Cahn-Hilliard equation.

ABSTRACT. Recent years have seen substantial interest in the development of high-order spatial discretizations for the Navier-Stokes equations, using either Scott-Vogelius elements (on suitable meshes) or H(div)-conforming elements to achieve high-order discretizations that strongly enforce the incompressibility constraint. For time-dependent problems, a further complication comes from achieving similar higher-order accuracy for the time-stepping scheme while maintaining optimal cost per time-step, roughly proportional to the number of spatial degrees of freedom. In this talk, we present simulation results using Firedrake for an H(div)-conforming spatial discretization and Irksome for fully implicit Runge-Kutta temporal discretizations, using PETSc’s PCPatch to define an effective block relaxation scheme for monolithic multigrid applied to the Jacobians of the resulting time-stepping scheme. We show that the scheme is robust in order, and discuss how the cost-per-timestep changes with spatial and temporal resolutions and discretization orders.

ABSTRACT. We investigate block preconditioners for large systems arising in implicit Runge-Kutta time stepping for hyperbolic partial differential equations. The system arising from IRK-time stepping applied to the wave equation is highly ill conditioned but has a block structure we exploit. If N is the number of degrees of freedom for the discretization of the PDE, then the resulting matrix after IRK-time stepping is a block s-by-s matrix where each block is of size 2N-by-2N. We present experimental results based on the LDU block preconditioner from Rana, Howle, Long, Meek, and Milestone (2021, SIAM J. Sci. Comp. v. 43, pp. 475-495), the block-Jacobi and the block Gauss-Siedel preconditioners from Mardal, Nilssen, and Staff (2007, SIAM J. Sci. Comp. v. 29, pp. 361-375), where we use GMRES for the stage equations and an AMG-V cycle for the subblocks. Our preliminary experiments show that LDU-based preconditioning is performing well in terms of reducing the condition number, improving the convergence rate and time as the spatial discretization is refined and as the number of RK stages increases.

ABSTRACT. WarpX (https://ecp-warpx.github.io/) is a massively parallel electromagnetic and electrostatic particle-in-cell code used for the modeling and simulation of a variety of plasma physics applications, including laser-plasma interactions, accelerator and beam physics, astrophysical plasmas, and microelectronics. WarpX scales to the world’s largest supercomputers and was awarded the 2022 ACM Gordon Bell Prize. In this talk we will present an overview of the WarpX code and its applications, with a special focus on the most relevant numerical methods and algorithms, including mesh refinement, load balancing, and multigrid for electrostatics. WarpX's underlying HPC software engineering strategy and roadmap will also be discussed.

ABSTRACT. It is well known that multigrid methods are very competitive in solving a wide range of SPD problems. However achieving such performance for non-SPD matrices remains an open problem. In particular, two main issues may arise when solving a Helmholtz problem. Some eigenvalues become negative or even complex, requiring the choice of an adapted smoothing method for capturing them. Moreover, since the near-kernel space is oscillatory, the geometric smoothness assumption cannot be used to build efficient interpolation rules. We present some investigations about designing a method that converges in a constant number of iterations with respect to the wavenumber. The method builds on an ideal reduction-based framework and related theory for SPD matrices to correct an initial least squares minimization coarse selection operator formed from a set of smoothed random vectors. We also present numerical results at the end of the paper.

ABSTRACT. We present a novel monolithic algebraic multigrid solver for the discrete Stokes problem discretized with stable mixed-finite elements. The algorithm is based on the use of a low-order $P_{1}isoP_{2}/P_{1}$ discretization as a preconditioner for a higher-order discretization, such as $P_2/P_1$. Smoothed aggregation algebraic multigrid is used to construct independent coarsening of the velocity and pressure fields, resulting in a purely algebraic preconditioner (i.e., using no geometric information). We show that this method effectively preconditions P_2/P_1$ (Taylor-Hood) and $P_2/P_1^{disc}$ (Scott-Vogelius) discretizations and can be further extended to higher-order discretizations.

ABSTRACT. Aggregate size and shape plays an important role in the overall convergence and complexity of smoothed-aggregation and root-node based algebraic multigrid methods. In this talk we present an algorithm that provides control over the number of aggregates and leads to ``well centered'' aggregates. Primarily we consider the challenge of forming balanced aggregates in the graph of a sparse matrix for use in multigrid, yet the algorithm has general applicability. Based on an extension of the Bellman-Ford algorithm, we generalize Lloyd's algorithm for partitioning a matrix graph to balance the number of nodes in each cluster. This is accompanied by a rebalancing algorithm that reduces the overall energy in the system. Theoretical results are provided to establish linear complexity and numerical results in the context of algebraic multigrid highlight the benefits of improved aggregation.

ABSTRACT. A crucial concept for algbraic multigrid (AMG) coarsening and interpolation is that of strength of connection (SoC) between degrees of freedom. The classical SoC measure is based on the relative sizes of entries in each row of the matrix and relies on heuristics that assume an M-matrix structure. This simple measure is cheap to evaluate and successful for many problems, but also relies on a user-defined strength threshold, which may need to be tuned for specific problems. In addition, the classical SoC measure can have difficulty identifying the proper strong connections for operators that are not M-matrices, particularly when there are both positive and negative off-diagonal entries in the same row. In this work, we examine low-cost techniques for building an auxiliary strength matrix that shares important properties with the original matrix operator while also being more amenable to the classical SoC measure. Applying classical SoC to this auxiliary strength matrix improves robustness of the SoC measure for a wider class of problems and across a wider range of strength thresholds.

ABSTRACT. It is well known that the classical algebraic multigrid (AMG) is a powerful iterative method for solving large linear systems that do not rely on the geometric structure of the problem. However, one of the limitations of AMG is that a pointwise smoother is typically assumed both theoretically and in practice. This assumption in turn limits the applicability of AMG leading to ineffective solvers for problems with a large near null-space. A canonical example is the Eddy-Current approximation of Maxwell Equations. To tackle this problem, block smoothers such as the distributive relaxation and the overlapping Schwarz smoothers are considered. In particular, we examine the effectiveness of compatible relaxation algorithms with these block smoothers in hope to develop a more robust general AMG algorithm.