Discretization Of Partial Differential Equations Research Articles

A randomized balanced proper orthogonal decomposition technique

In this paper, we consider the model reduction problem of large-scale linear systems, such as systems obtained through the discretization of partial differential equations. We propose a randomized proper orthogonal decomposition (RPOD∗) technique to obtain the reduced order model by perturbing the primal and adjoint system using Gaussian white noise. We show that computations required by RPOD∗ algorithm is orders of magnitude lower while its performance is much better than other state of the art algorithms. We also relate the RPOD∗ algorithm to Krylov subspace methods and show that it constitutes a randomized approach to computational linear algebra problems that utilize Krylov subspace methods.

Extension of Tensor-Product Generalized and Dense-Norm Summation-by-Parts Operators to Curvilinear Coordinates

Methodologies are presented that enable the construction of provably linearly stable and conservative high-order discretizations of partial differential equations in curvilinear coordinates based on generalized summation-by-parts operators, including operators with dense-norm matrices. Specifically, three approaches are presented for the construction of stable and conservative schemes in curvilinear coordinates using summation-by-parts (SBP) operators that have a diagonal norm but may or may not include boundary nodes: (1) the mortar-element approach, (2) the global SBP-operator approach, and (3) the staggered-grid approach. Moreover, the staggered-grid approach is extended to enable the development of stable dense-norm operators in curvilinear coordinates. In addition, collocated upwind simultaneous approximation terms for the weak imposition of boundary conditions or inter-element coupling are extended to curvilinear coordinates with the new approaches. While the emphasis in the paper is on tensor-product SBP operators, the approaches that are covered are directly applicable to multidimensional SBP operators.

Learning data-driven discretizations for partial differential equations

The numerical solution of partial differential equations (PDEs) is challenging because of the need to resolve spatiotemporal features over wide length- and timescales. Often, it is computationally intractable to resolve the finest features in the solution. The only recourse is to use approximate coarse-grained representations, which aim to accurately represent long-wavelength dynamics while properly accounting for unresolved small-scale physics. Deriving such coarse-grained equations is notoriously difficult and often ad hoc. Here we introduce data-driven discretization, a method for learning optimized approximations to PDEs based on actual solutions to the known underlying equations. Our approach uses neural networks to estimate spatial derivatives, which are optimized end to end to best satisfy the equations on a low-resolution grid. The resulting numerical methods are remarkably accurate, allowing us to integrate in time a collection of nonlinear equations in 1 spatial dimension at resolutions 4× to 8× coarser than is possible with standard finite-difference methods.

Radial discrete PDE splines on Lipschitz domains

In the framework of an elliptic partial differential equation (PDE), certain boundary conditions and a set of points to approximate in a Lipschitz domain and arbitrary dimension, we use radial basis function (RBF) techniques for the construction and characterization of discrete PDE splines. We also show convergence and derive error estimates.

The INTERNODES Method for Non-conforming Discretizations of PDEs

INTERNODES is a general purpose method to deal with non-conforming discretizations of partial differential equations on 2D and 3D regions partitioned into two or several disjoint subdomains. It exploits two intergrid interpolation operators, one for transfering the Dirichlet trace across the interfaces, and the other for the Neumann trace. In this paper, in every subdomain the original problem is discretized by either the finite element method (FEM) or the spectral element method (SEM or hp-FEM), using a priori non-matching grids and piecewise polynomials of different degrees. Other discretization methods, however, can be used. INTERNODES can also be applied to heterogeneous or multiphysics problems, that is, problems that feature different differential operators inside adjacent subdomains. For instance, in this paper we apply the INTERNODES method to a Stokes–Darcy coupled problem that models the filtration of fluids in porous media. Our results highlight the flexibility of the method as well as its optimal rate of convergence with respect to the grid size and the polynomial degree.

Novel multilevel techniques for convergence acceleration in the solution of systems of equations arising from RBF-FD meshless discretizations

The present paper develops two new techniques, namely additive correction multicloud (ACMC) and smoothed restriction multicloud (SRMC), for the efficient solution of systems of equations arising from Radial Basis Function-generated Finite Difference (RBF-FD) meshless discretizations of partial differential equations (PDEs). RBF-FD meshless methods employ arbitrary distributed nodes, without the need to generate a mesh, for the numerical solution of PDEs. The proposed techniques are specifically designed for the RBF-FD data structure and employ simple restriction and interpolation strategies in order to obtain a hierarchy of coarse-level node distributions and the corresponding correction equations.Both techniques are kept as simple as possible in terms of code implementation, which is an important feature of meshless methods. The techniques are verified on 2D and 3D Poisson equations, defined on non-trivial domains , showing very high benefits in terms of both time consumption and work to convergence when comparing the present techniques to the most common solver approaches. These benefits make the RBF-FD approach competitive with standard grid-based approaches when the number of nodes is very high, allowing large size problems to be tackled by the RBF-FD method.

Robust approaches to handling complex geometries with Galerkin difference methods

The Galerkin difference (GD) basis is a set of continuous, piecewise polynomials defined using a finite-difference-like grid of degrees of freedom. The one dimensional GD basis functions are naturally extended to multiple dimensions using a tensor product construction on quadrilateral elements. The GD basis can be used to define the solution space for a discontinuous Galerkin finite element discretization of partial differential equations. In this work we propose two approaches to handling complex geometries within this setting: (1) using nonconforming, curvilinear GD elements and (2) coupling affine GD elements with curvilinear simplicial elements. In both cases the (semidiscrete) discontinuous Galerkin method is provably energy stable even when variational crimes are committed. Additionally, for both element types a weight-adjusted mass matrix is used, which ensures that only the reference mass matrix must be inverted. We also present sufficient conditions on the treatment of metric terms for the curvilinear, nonconforming GD elements to ensure that the scheme is both constant preserving and conservative. Numerical experiments confirm the stability results and demonstrate the accuracy of the coupled schemes.

AMGCL: An Efficient, Flexible, and Extensible Algebraic Multigrid Implementation

The paper presents AMGCL -- an opensource C++ library implementing the algebraic multigrid method (AMG) for solution of large sparse linear systems of equations, usually arising from discretization of partial differential equations on an unstructured grid. The library supports both shared and distributed memory computation, allows to utilize modern massively parallel processors via OpenMP, OpenCL, or CUDA technologies, has minimal dependencies, and is easily extensible. The design principles behind AMGCL are discussed and it is shown that the code performance is on par with alternative implementations.

A localized reduced-order modeling approach for PDEs with bifurcating solutions

Reduced-order modeling (ROM) commonly refers to the construction, based on a few solutions (referred to as snapshots) of an expensive discretized partial differential equation (PDE), and the subsequent application of low-dimensional discretizations of partial differential equations (PDEs) that can be used to more efficiently treat problems in control and optimization, uncertainty quantification, and other settings that require multiple approximate PDE solutions. Although ROMs have been successfully used in many settings, ROMs built specifically for the efficient treatment of PDEs having solutions that bifurcate as the values of input parameters change have not received much attention. In such cases, the parameter domain can be subdivided into subregions, each of which corresponds to a different branch of solutions. Popular ROM approaches such as proper orthogonal decomposition (POD), results in a global low-dimensional basis that does not respect the often large differences in the PDE solutions corresponding to different subregions. In this work, we develop and test a new ROM approach specifically aimed at bifurcation problems. In the new method, the k-means algorithm is used to cluster snapshots so that within cluster snapshots are similar to each other and are dissimilar to those in other clusters. This is followed by the construction of local POD bases, one for each cluster. The method also can detect which cluster a new parameter point belongs to, after which the local basis corresponding to that cluster is used to determine a ROM approximation. Numerical experiments show the effectiveness of the method both for problems for which bifurcation cause continuous and discontinuous changes in the solution of the PDE.

Approximation properties of hybrid shearlet-wavelet frames for Sobolev spaces

In this paper, we study a newly developed shearlet system on bounded domains which yields frames for $H^s(\Omega)$ for some $s\in \mathbb{N}$, $\Omega \subset \mathbb{R}^2$. We will derive approximation rates with respect to $H^s(\Omega)$ norms for functions whose derivatives admit smooth jumps along curves and demonstrate superior rates to those provided by pure wavelet systems. These improved approximation rates demonstrate the potential of the novel shearlet system for the discretization of partial differential equations. Therefore, we implement an adaptive shearlet-based algorithm for the solution of an elliptic PDE and analyze its computational complexity and convergence properties.

Nonlinear Circle-Criterion Observer Design for an Electrochemical Battery Model

Smart battery management systems require reliable information on the internal state of the battery, which are unavailable through direct measurements. Electrochemical models are relevant in this context as these describe the internal phenomena governing the battery dynamics. These models can thus be used to design observers, which provide estimates of the state variables online. In this paper, we consider a single particle electrochemical model of a lithium-ion battery given by a set of ordinary differential equations built from the spatial discretization of partial differential equations (PDEs) that locally describe the mass and charge transport of lithium. We then design a nonlinear observer, whose global exponential stability and robustness are guaranteed provided that a linear matrix inequality holds. The latter is shown to be verified for standard model parameters values. Simulation results based on plug-in hybrid electric vehicle tests on the original PDE model are presented, which illustrate the good performances of the observer in terms of speed of convergence and robustness.

Positivity Preserving Limiters for Time-Implicit Higher Order Accurate Discontinuous Galerkin Discretizations

Currently, nearly all positivity preserving discontinuous Galerkin (DG) discretizations of partial differential equations are coupled with explicit time integration methods. Unfortunately, for many problems this can result in severe time-step restrictions. The techniques used to develop explicit positivity preserving DG discretizations cannot, however, easily be combined with implicit time integration methods. In this paper, we therefore present a new approach. Using Lagrange multipliers, the conditions imposed by the positivity preserving limiters are directly coupled to a DG discretization combined with a diagonally implicit Runge--Kutta time integration method. The positivity preserving DG discretization is then reformulated as a Karush--Kuhn--Tucker (KKT) problem, which is frequently encountered in constrained optimization. Since the limiter is only active in areas where positivity must be enforced, it does not affect the higher order DG discretization elsewhere. The resulting nonsmooth nonlinear algebraic equations have, however, a different structure compared to most constrained optimization problems. We therefore develop an efficient active set semismooth Newton method that is suitable for the KKT formulation of time-implicit positivity preserving DG discretizations. Convergence of this semismooth Newton method is proven using a specially designed quasi-directional derivative of the time-implicit positivity preserving DG discretization. The time-implicit positivity preserving DG discretization is demonstrated for several nonlinear scalar conservation laws, which include the advection, Burgers, Allen--Cahn, Barenblatt, and Buckley--Leverett equations.

Sparse supernodal solver using block low-rank compression: Design, performance and analysis

This paper presents two approaches using a Block Low-Rank (BLR) compression technique to reduce the memory footprint and/or the time-to-solution of the sparse supernodal solver PaStiX. This flat, non-hierarchical, compression method allows to take advantage of the low-rank property of the blocks appearing during the factorization of sparse linear systems, which come from the discretization of partial differential equations. The proposed solver can be used either as a direct solver at a lower precision or as a very robust preconditioner. The first approach, called Minimal Memory, illustrates the maximum memory gain that can be obtained with the BLR compression method, while the second approach, called Just-In-Time, mainly focuses on reducing the computational complexity and thus the time-to-solution. Singular Value Decomposition (SVD) and Rank-Revealing QR (RRQR), as compression kernels, are both compared in terms of factorization time, memory consumption, as well as numerical properties. Experiments on a shared memory node with 24 threads and 128 GB of memory are performed to evaluate the potential of both strategies. On a set of matrices from real-life problems, we demonstrate a memory footprint reduction of up to 4 times using the Minimal Memory strategy and a computational time speedup of up to 3.5 times with the Just-In-Time strategy. Then, we study the impact of configuration parameters of the BLR solver that allowed us to solve a 3D laplacian of 36 million unknowns a single node, while the full-rank solver stopped at 8 million due to memory limitation.

Denoising of Image Gradients and Total Generalized Variation Denoising

We revisit total variation denoising and study an augmented model where we assume that an estimate of the image gradient is available. We show that this increases the image reconstruction quality and derive that the resulting model resembles the total generalized variation denoising method, thus providing a new motivation for this model. Further, we propose to use a constraint denoising model and develop a variational denoising model that is basically parameter free, i.e., all model parameters are estimated directly from the noisy image. Moreover, we use Chambolle–Pock’s primal dual method as well as the Douglas–Rachford method for the new models. For the latter one has to solve large discretizations of partial differential equations. We propose to do this in an inexact manner using the preconditioned conjugate gradients method and derive preconditioners for this. Numerical experiments show that the resulting method has good denoising properties and also that preconditioning does increase convergence speed significantly. Finally, we analyze the duality gap of different formulations of the TGV denoising problem and derive a simple stopping criterion.

Efficient Monte Carlo With Graph‐Based Subsurface Flow and Transport Models

AbstractSimulating flow and transport in fractured porous media frequently involves solving numerical discretizations of partial differential equations with a large number of degrees of freedom using discrete fracture network (DFN) models. Uncertainty in the properties of the fracture network that controls flow and transport requires a large number of DFN simulations to statistically describe quantities of interest. However, the computational cost of solving more than a few realizations of a large DFN can be intractable. As a means of circumventing this problem, we utilize both a high‐fidelity DFN model and a graph‐based model of flow and transport in combination with a multifidelity Monte Carlo (MC) method to reduce the number of high‐fidelity simulations that are needed to obtain an accurate estimate of the quantity of interest. We demonstrate the approach by estimating quantiles of the breakthrough time for a conservative tracer in an ensemble of fractured porous media. Our results demonstrate that a multifidelity MC estimate, whose computational cost is equal to the cost of 10 DFN simulations, can be as accurate as a standard MC estimate that utilizes 1,000 DFN simulations. Thus the combination of our graph‐based model with multifidelity MC estimates effectively reduces the computational cost of the problem by a factor of approximately 100.

Partial tensor decomposition for decoupling isogeometric Galerkin discretizations

System matrix assembly for isogeometric (i.e., spline-based) discretizations of partial differential equations is more challenging than for classical finite elements, due to the increased polynomial degrees and the larger (and hence more overlapping) supports of the basis functions. The global tensor-product structure of the discrete spaces employed in isogeometric analysis can be exploited to accelerate the computations, using sum factorization, precomputed look-up tables, and tensor decomposition. We generalize the third approach by considering partial tensor decompositions. We show that the resulting new method preserves the global discretization error and that its computational complexity compares favorably to the existing approaches. Moreover, the numerical realization simplifies considerably since it relies on standard techniques from numerical linear algebra.

Numerical low-rank approximation of matrix differential equations

Abstract The efficient numerical integration of large-scale matrix differential equations is a topical problem in numerical analysis and of great importance in many applications. Standard numerical methods applied to such problems require an unduly amount of computing time and memory, in general. Based on a dynamical low-rank approximation of the solution, a new splitting integrator is proposed for a quite general class of stiff matrix differential equations. This class comprises differential Lyapunov and differential Riccati equations that arise from spatial discretizations of partial differential equations. The proposed integrator handles stiffness in an efficient way, and it preserves the symmetry and positive semidefiniteness of solutions of differential Lyapunov equations. Numerical examples that illustrate the benefits of this new method are given. In particular, numerical results for the efficient simulation of the weather phenomenon El Nino are presented.

Sparsity preserving optimal control of discretized PDE systems

We focus on the problem of optimal control of large-scale systems whose models are obtained by discretization of partial differential equations using the Finite Element (FE) or Finite Difference (FD) methods. The motivation for studying this pressing problem originates from the fact that the classical numerical tools used to solve low-dimensional optimal control problems are computationally infeasible for large-scale systems. Furthermore, although the matrices of large-scale FE or FD models are usually sparse banded or highly structured, the optimal control solution computed using the classical methods is dense and unstructured. Consequently, it is not suitable for efficient centralized and distributed real-time implementations. We show that the a priori (sparsity) patterns of the exact solutions of the generalized Lyapunov equations for FE and FD models are banded matrices. The a priori pattern predicts the dominant non-zero entries of the exact solution. We furthermore show that for well-conditioned problems, the a priori patterns are not only banded but also sparse matrices. On the basis of these results, we develop two computationally efficient methods for computing sparse approximate solutions of generalized Lyapunov equations. Using these two methods and the inexact Newton method, we show that the solution of the generalized Riccati equation can be approximated by a banded matrix. This enables us to develop a novel computationally efficient optimal control approach that is able to preserve the sparsity of the control law. We perform extensive numerical experiments that demonstrate the effectiveness of our approach.

An Enriched Multiscale Mortar Space for High Contrast Flow Problems

Mortar methods are widely used techniques for discretizations of partial differential equations and preconditioners for the algebraic systems resulting from the discretizations. For problems with high contrast and multiple scales, the standard mortar spaces are not robust, and some enrichments are necessary in order to obtain an efficient and robust mortar space. In this paper, we consider a class of flow problems in high contrast heterogeneous media, and develop a systematic approach to obtain an enriched multiscale mortar space. Our approach is based on the constructions of local multiscale basis functions. The multiscale basis functions are constructed from local problems by following the framework of the Generalized Multiscale Finite Element Method (GMsFEM). In particular, we first create a local snapshot space. Then we select the dominated modes within the snapshot space using an appropriate Proper Orthogonal Decomposition (POD) technique. These multiscale basis functions show better accuracy than polynomial basis for multiscale problems. Using the proposed multiscale mortar space, we will construct a multiscale finite element method to solve the flow problem on a coarse grid and a preconditioning technique for the fine scale discretization of the flow problem. In particular, we develop a multiscale mortar mixed finite element method using the mortar space. In addition, we will design a two-level additive preconditioner and a two-level hybrid preconditioner based on the proposed mortar space for the iterative method applied to the fine scale discretization of the flow problem. We present several numerical examples to demonstrate the efficiency and robustness of our proposed mortar space with respect to both the coarse multiscale solver and the preconditioners.

Asynchronous Parareal Time Discretization For Partial Differential Equations

Asynchronous iterations have been investigated more and more for both scaling and fault-resilience purposes on high performance computing platforms. While so far, they have been exclusively applied within space domain decomposition frameworks, this paper advocates a novel application direction targeting time-decomposed time-parallel approaches. Specifically, an asynchronous iterative model is derived from the Parareal scheme, for which convergence and speedup analysis are then conducted. It turned out that Parareal and async-Parareal feature very close convergence conditions, asymptotically equivalent, including the finite-time termination property. Based on a computational cost model aware of unsteady communication delays, our speedup analysis shows the potential performance gain from asynchronous iterations, which is confirmed by some experimental cases of heat evolution on a homogeneous supercomputer. This primary work clearly suggests possible further benefits from asynchronous iterations.