Consider an algebraic two-level method applied to the n-dimensional linear system A\bfitx = \bfitbusing fine-space preconditioner (i.e., "relaxation" or "smoother") M, with M \approxA, restriction and interpolation R and P, and algebraic coarse-space operator A(c) := R*AP. Then, what are the best possible transfer operators R and P of a given dimension nc < n? Brannick et al. [12] showed that when A and M are Hermitian positive definite (HPD), the optimal interpolation is such that its range contains the nc smallest generalized eigenvectors of the matrix pencil (A, M). Recently, in Ali et al. [5] we generalized this framework to the non-HPD setting, by considering both right (interpolation) and left (restriction) generalized eigenvectors of (A, M) and defining corresponding nonsymmetric transfer operators {R-#, P-#}. Tight convergence bounds for {R-#, P-#} are derived in spectral radius, as well as a proof of pseudo-optimality. Note, {R-#, P-#} are typically complex valued, which is not practical for real-valued problems. Here we build on [5], first characterizing all inner products in which the coarse-space correction defined by {R-#, P-#} is orthogonal. We then develop tight two-level convergence bounds in these norms, and prove that the underlying transfer operators {R-#, P-#} are genuinely optimal. As a special case, our theory both recovers and extends the HPD results from [12]. Finally, we show how to construct optimal, real-valued transfer operators in the case of that A and M are real valued, but are not HPD. Numerical examples arising from a discretized advection-reaction equation, wave-equation, and Stokes equations are used to verify and illustrate the theory.

We propose a sparse interpolation construction and a practical coarsening algorithm for the algebraic multigrid (AMG) method, tailored towards H(curl). Building on the generalized AMG framework, we introduce an interior/exterior splitting that yields both a refinement-based and a fully algebraic construction of the interpolation. The refinement-based approach follows geometric hierarchy, while the purely algebraic interpolation is constructed through a coarsening process that first coarsens a nodal dual problem and then builds coarse and fine variables using a matching algorithm. We establish the weak approximation property and the commuting relation under certain assumptions. Combined with matching block smoothers, the proposed interpolation yields an effective algebraic multilevel method. Numerical experiments show robustness under strong coefficient jumps, where the proposed methods substantially outperform standard geometric multigrid.

We present an approach to constructing a practical coarsening algorithm and interpolation operator for the algebraic multigrid (AMG) method, tailored towards systems of partial differential equations (PDEs) with large near-kernels, such as H(curl) and H(div). Our method builds on compatible relaxation (CR) and the ideal interpolation model within the generalized AMG (GAMG) framework but introduces several modifications to define an AMG method for PDE systems. We construct an interpolation operator through a coarsening process that first coarsens a nodal dual problem and then builds the coarse and fine variables using a matching algorithm. Our interpolation follows the ideal formulation; however, we enhance the sparsity of ideal interpolation by decoupling the fine and coarse variables completely. When the coarse variables align with the geometric refinement, our method reproduces re-discretization on unstructured meshes. Together with an automatic smoother construction scheme that identifies the local near kernels, our approach forms a complete two-grid method. Finally, we also show numerical results that demonstrate the effectiveness of this interpolation scheme by applying it to targeted problems and the Stokes system.

We present multiscale graph-based reduction algorithms for upscaling heterogeneous and anisotropic diffusion problems. The proposed coarsening approaches begin by constructing a partitioning of the computational domain into a set of balanced local subdomains, resulting in a standard type of domain decomposition. Given this initial decomposition, general coarsening techniques based on spectral clustering are applied within each subgraph in order to accurately identify the key microscopic features of a given system. The spectral clustering algorithm is based on local generalized eigen-decompositions applied to the signed graph Laplacian. The resulting coarse-fine splittings are combined with two variants of energy-minimizing strategies for constructing coarse bases for diffusion problems. The first is an unconstrained minimization formulation in which local harmonic extensions are applied column-wise to construct multi-vector preserving interpolation in each region, whereas the second approach is a variant of the constrained energy minimization formulations derived in the context of non-local multi-continua upscaling techniques. We apply the resulting upscaling algorithms to a variety of tests coming from the graph Laplacian, including diffusion in the perforated domain, channelized media, highly anisotropic settings, and discrete pore network models to demonstrate the potential and robustness of the proposed coarsening approaches. We show numerically and theoretically that the proposed approaches lead to accurate coarse-scale models.

Algebraic multigrid (AMG) is known to be an effective solver for many sparse symmetric positive definite (SPD) linear systems. For SPD systems, the convergence theory of AMG is well-understood in terms of the A-norm, but in a nonsymmetric setting, such an energy norm is non-existent. For this reason, convergence of AMG for nonsymmetric systems of equations remains an open area of research. A particular aspect missing from theory of nonsymmetric and indefinite AMG is the incorporation of general relaxation schemes. In the SPD setting, the classical form of optimal AMG interpolation provides a useful insight in determining the best possible two-grid convergence rate of a method based on an arbitrary symmetrized relaxation scheme. In this work, we discuss a generalization of the optimal AMG convergence theory targeting nonsymmetric problems, using a certain matrix-induced orthogonality of the left and right eigenvectors of a generalized eigenvalue problem relating the system matrix and relaxation operator. We show that using this generalization of the optimal convergence theory, one can obtain a measure of the spectral radius of the two grid error transfer operator that is mathematically equivalent to the derivation in the SPD setting for optimal interpolation, which instead uses norms. In addition, this generalization of the optimal AMG convergence theory can be further extended for symmetric indefinite problems, such as those arising from saddle point systems so that one can obtain a precise convergence rate of the resulting two-grid method based on optimal interpolation. We provide supporting numerical examples of the convergence theory for nonsymmetric advection-diffusion problems, two-dimensional Dirac equation motivated by γ_5-symmetry, and the mixed Darcy flow problem corresponding to a saddle point system.

This paper focuses on developing a reduction-based algebraic multigrid method that is suitable for solving general (non)symmetric linear systems and is naturally robust from pure advection to pure diffusion. Initial motivation comes from a new reduction-based algebraic multigrid (AMG) approach, $\ell$AIR (local approximate ideal restriction), that was developed for solving advection-dominated problems. Though this new solver is very effective in the advection dominated regime, its performance degrades in cases where diffusion becomes dominant. This is consistent with the fact that in general, reduction-based AMG methods tend to suffer from growth in complexity and/or convergence rates as the problem size is increased, especially for diffusion dominated problems in two or three dimensions. Motivated by the success of $\ell$AIR in the advective regime, our aim in this paper is to generalize the AIR framework with the goal of improving the performance of the solver in diffusion dominated regimes. To do so, we propose a novel way to combine mode constraints as used commonly in energy minimization AMG methods with the local approximation of ideal operators used in $\ell$AIR. The resulting constrained $\ell$AIR (C$\ell$AIR) algorithm is able to achieve fast scalable convergence on advective and diffusive problems. In addition, it is able to achieve standard low complexity hierarchies in the diffusive regime through aggressive coarsening, something that has been previously difficult for reduction-based methods.

The Pennsylvania State University (“Subcontractor”) worked on developing new parallel algebraic multilevel methods suitable for solving PDEs. Specifically, work on the design of multigrid solvers for coupled systems of partial differential equations arising in numerical modeling of various applications was completed. A main emphasis was on the design of new ideal algebraic multigrid interpolation for problems such as Maxwell’s equations where block smoothers are needed and the standard form of ideal interpolation is not an effective choice.

This paper introduces bootstrap multigrid methods for solving eigenvalue problems arising from the discretization of partial differential equations. Inspired by the full bootstrap algebraic multigrid (BAMG) setup algorithm that includes an AMG eigensolver, it is illustrated how the algorithm can be simplified for the case of a discretized partial differential equation (PDE), thereby developing a bootstrap geometric multigrid (BMG) approach. We illustrate numerically the efficacy of the BMG method for: (1) recovering eigenvalues having large multiplicity, (2) computing interior eigenvalues, and (3) approximating shifted indefinite eigenvalue problems. Numerical experiments are presented to illustrate the basic components and ideas behind the success of the overall bootstrap multigrid approach. For completeness, we present a simplified error analysis of a two-grid bootstrap algorithm for the Laplace-Beltrami eigenvalue problem.

The project team continued with work on developing new bootstrap AMG techniques for solving symmetric and non-symmetric PDE systems. The focus of this work is to derive more reliable measures of the quality of the coarse space set than the convergence rate of the standard F-relaxation form of CR and a more robust form of interpolation than the so-called ideal form. We have successfully derived a sharp variant of CR that gives the precise convergence rate of the two-grid method using this optimal interpolation and, in addition, we derived a new Generalized Bootstrap AMG setup algorithm that uses as its main tool a multilevel eigensolver for the generalized eigenvalue problem involving the system matrix and the symmetrized smoother. In addition, the approach allows for general block smoothers with overlap. We have applied the method to scalar diffusion problems, linear elasticity, and Maxwell’s and the method shows marked improvements over existing AMG methods for these problems. In addition, the team worked with CASC members on new forms of ideal AMG interpolation.

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

In this paper, two approaches for modeling three-component fluid flows using diffusive interface method are discussed. Thermodynamic consistency of the proposed models is preserved when using an energetic variational framework to derive the coupled systems of partial differential equations that comprise the resulting models. The issue of algebraic and dynamic consistency is investigated. In addition, the two approaches that are presented are compared analytically and numerically.

Algebraic multigrid (AMG) methods are powerful solvers with linear or near-linear computational complexity for certain classes of linear systems, Ax=b. Broadening the scope of problems that AMG can effectively solve requires the development of improved interpolation operators. Such development is often based on AMG convergence theory. However, convergence theory in AMG tends to have a disconnect with AMG in practice due to the practical constraints of (i) maintaining matrix sparsity in transfer and coarse-grid operators, and (ii) retaining linear complexity in the setup and solve phase. This paper presents a review of fundamental results in AMG convergence theory, followed by a discussion on how these results can be used to motivate interpolation operators in practice. A general weighted energy minimization functional is then proposed to form interpolation operators, and a novel `diagonal' preconditioner for Sylvester- or Lyapunov-type equations developed simultaneously. Although results based on the weighted energy minimization typically underperform compared to a fully constrained energy minimization, numerical results provide new insight into the role of energy minimization and constraint vectors in AMG interpolation.

Progress has been made in several of the proposed areas of research, as summarized below: Compatible Relaxation and Optimal Interpolation and Bootstrap Algebraic Multigrid and non-Galerkin AMG for the Dirac equations. Our preliminary findings are promising and suggest that these techniques should lead to efficient and robust solvers for various discretizations of the Dirac PDE.

In this paper, we consider a classical algebraic multigrid (AMG) form of optimal interpolation that directly minimizes the two-grid convergence rate and compare it with a so-called ideal interpolation that minimizes a weak approximation property of the coarse space. We study compatible relaxation type estimates for the quality of the coarse grid and derive a new sharp measure using optimal interpolation that provides a guaranteed lower bound on the convergence rate of the resulting two-grid method for a given grid. In addition, we design a generalized bootstrap AMG setup algorithm that computes a sparse approximation to the optimal interpolation matrix. We demonstrate numerically that the bootstrap AMG method with sparse interpolation matrix (and spanning multiple levels) converges faster than the two-grid method with the standard ideal interpolation (a dense matrix) for various scalar diffusion problems with highly varying diffusion coefficient.

to construct geologically constrained models of faulted/fractured reservoirs. Discrete Fracture Network (DFN) simulation provides the option of performing multiphase calculations on spatially explicit, geologically feasible fracture sets. Multiphase DFN simulations of and sensitivity studies on a wide variety of fracture networks created using fracture creation/simulation programs was undertaken in the first part of this project. This involved creating interfaces to seamlessly convert the fracture characterization information into simulator input, grid the complex geometry, perform the simulations, and analyze and visualize results. Benchmarking and comparison with conventional simulators was also a component of this work. After demonstration of the fact that multiphase simulations can be carried out on complex fracture networks, quantitative effects of the heterogeneity of fracture properties were evaluated. Reservoirs are populated with fractures of several different scales and properties. A multiscale fracture modeling study was undertaken and the effects of heterogeneity and storage on water displacement dynamics in fractured basements were investigated. In gravity-dominated systems, more oil could be recovered at a given pore volume of injection at lower rates. However, if oil production can be continued at high water cuts, the discounted cumulative production usually favors higher production rates. The workflow developed during the project was also used to perform multiphase simulations in heterogeneous, fracture-matrix systems. Compositional and thermal-compositional simulators were developed for fractured reservoirs using the generalized framework. The thermal-compositional simulator was based on a novel 'equation-alignment' approach that helped choose the correct variables to solve depending on the number of phases present and the prescribed component partitioning. The simulators were used in steamflooding and in insitu combustion applications. The framework was constructed to be inherently parallel. The partitioning routines employed in the framework allowed generalized partitioning on highly complex fractured reservoirs and in instances when wells (incorporated in these models as line sources) were divided between two or more processors.

The diffusive interface method is an approach for modeling interactions among complex substances. The main idea behind this method is to introduce phase field labeling functions in order to model the contact line by smooth change from one type of material to another. The approach has been widely used and successfully incorporated to numerous practical applications, including models of phase transitions, contact line dynamics in complex fluids, cell motility. and many other problems in science and engineering. In this review article, we describe how the diffusive interface method can be used for the modeling of complex fluids and point to some main results in the field.

The overlap operator is a lattice discretization of the Dirac operator of quantum chromodynamics (QCD), the fundamental physical theory of the strong interaction between the quarks. As opposed to other discretizations, it preserves the important physical property of chiral symmetry, at the expense of requiring much more effort when solving systems posed with this operator. We present a preconditioning technique based on another lattice discretization, the Wilson-Dirac operator. The mathematical analysis precisely describes the effect of this preconditioning strategy in the case that the Wilson-Dirac operator is normal. Although this is not exactly the case in realistic settings, we show that current smearing techniques indeed drive the Wilson-Dirac operator towards normality, thus providing motivation for why our preconditioner works well in practice. Results of numerical experiments in physically relevant settings show that our preconditioning yields accelerations of more than an order of magnitude compared to unpreconditioned solvers.

The project team has continued with work on developing aggressive coarsening techniques for AMG methods. Of particular interest is the idea to use aggressive coarsening with polynomial smoothing. Using local Fourier analysis the optimal values for the parameters involved in defining the polynomial smoothers are determined automatically in a way to achieve fast convergence of cycles with aggressive coarsening. Numerical tests have the sharpness of the theoretical results. The methods are highly parallelizable and efficient multigrid algorithms on structured and semistructured grids in two and three spatial dimensions.