Showing papers in "SIAM Journal on Scientific Computing in 2018"
TL;DR: In this paper, a new algorithm for approximation by rational functions on a real or complex set of points, implementable in 40 lines of MATLAB and requiring no user input parameters, is presented.
Abstract: We introduce a new algorithm for approximation by rational functions on a real or complex set of points, implementable in 40 lines of MATLAB and requiring no user input parameters. Even on a disk o...
271 citations
TL;DR: The method circumvents the need for spatial discretization of the differential operators by proper placement of Gaussian process priors and is an attempt to construct structured and data-efficient learning machines, which are explicitly informed by the underlying physics that possibly generated the observed data.
Abstract: We introduce the concept of numerical Gaussian processes, which we define as Gaussian processes with covariance functions resulting from temporal discretization of time-dependent partial differential equations. Numerical Gaussian processes, by construction, are designed to deal with cases where (a) all we observe are noisy data on black-box initial conditions, and (b) we are interested in quantifying the uncertainty associated with such noisy data in our solutions to time-dependent partial differential equations. Our method circumvents the need for spatial discretization of the differential operators by proper placement of Gaussian process priors. This is an attempt to construct structured and data-efficient learning machines, which are explicitly informed by the underlying physics that possibly generated the observed data. The effectiveness of the proposed approach is demonstrated through several benchmark problems involving linear and nonlinear time-dependent operators. In all examples, we are able to r...
218 citations
TL;DR: For solving large-scale systems of linear equations by iteration methods, an effective probability criterion for selecting the working rows from the coefficient matrix is introduced.
Abstract: For solving large-scale systems of linear equations by iteration methods, we introduce an effective probability criterion for selecting the working rows from the coefficient matrix and construct a ...
154 citations
TL;DR: In this article, a shift proper orthogonal decomposition (sPOD) is proposed for transport-dominated model reduction, which features a determination of the multiple transport velocities and a separation of these.
Abstract: Transport-dominated phenomena provide a challenge for common mode-based model reduction approaches. We present a model reduction method, which is suited for these kinds of systems. It extends the proper orthogonal decomposition (POD) by introducing time-dependent shifts of the snapshot matrix. The approach, called shifted proper orthogonal decomposition (sPOD), features a determination of the multiple transport velocities and a separation of these. One- and two-dimensional test examples reveal the good performance of the sPOD for transport-dominated phenomena and its superiority in comparison to the POD.
145 citations
TL;DR: The results suggest that on architectures for which half precision is efficiently implemented it will be possible to solve certain linear systems up to twice as fast and to greater accuracy, as well as recommending a standard solver that uses LU factorization in single precision.
Abstract: We propose a general algorithm for solving a $n\times n$ nonsingular linear system $Ax = b$ based on iterative refinement with three precisions. The working precision is combined with possibly different precisions for solving for the correction term and for computing the residuals. Via rounding error analysis of the algorithm we derive sufficient conditions for convergence and bounds for the attainable normwise forward error and normwise and componentwise backward errors. Our results generalize and unify many existing rounding error analyses for iterative refinement. With single precision as the working precision, we show that by using LU factorization in IEEE half precision as the solver and calculating the residuals in double precision it is possible to solve $Ax = b$ to full single precision accuracy for condition numbers $\kappa_2(A) \le 10^4$, with the $O(n^3)$ part of the computations carried out entirely in half precision. We show further that by solving the correction equations by GMRES preconditioned by the LU factors the restriction on the condition number can be weakened to $\kappa_2(A) \le 10^8$, although in general there is no guarantee that GMRES will converge quickly. Taking for comparison a standard $Ax = b$ solver that uses LU factorization in single precision, these results suggest that on architectures for which half precision is efficiently implemented it will be possible to solve certain linear systems $Ax = b$ up to twice as fast \emph{and} to greater accuracy. Analogous results are given with double precision as the working precision.
144 citations
TL;DR: In this article, a data-driven filtered reduced order model (DDF-ROM) framework is proposed for numerical simulation of fluid flows, which is based on general ideas of spatial filtering and optimization and is independent of (restrictive) phenomenological arguments.
Abstract: We propose a data-driven filtered reduced order model (DDF-ROM) framework for the numerical simulation of fluid flows. The novel DDF-ROM framework consists of two steps: (i) In the first step, we use ROM projection to filter the nonlinear PDE and construct a filtered ROM. This filtered ROM is low-dimensional but is not closed (because of the nonlinearity in the given PDE). (ii) In the second step, we use data-driven modeling to close the filtered ROM, i.e., to model the interaction between the resolved and unresolved modes. To this end, we use a quadratic ansatz to model this interaction and close the filtered ROM. To find the new coefficients in the closed filtered ROM, we solve an optimization problem that minimizes the difference between the full order model data and our ansatz. We emphasize that the new DDF-ROM is built on general ideas of spatial filtering and optimization and is independent of (restrictive) phenomenological arguments. We investigate the DDF-ROM in the numerical simulation of a 2D ch...
129 citations
TL;DR: The multiple scalar auxiliary variable (MSAV) approach is applied to the phase-field vesicle membrane (PF-VMEM) model which, in addition to some usual nonlinear terms in the free energy, has two additional penalty terms to enforce the volume and surface area.
Abstract: We consider in this paper gradient flows with disparate terms in the free energy that cannot be efficiently handled with the scalar auxiliary variable (SAV) approach, and we develop the multiple sc...
82 citations
TL;DR: The linearization technique developed in this paper is so general that it can be applied to any thermodynamically consistent hydrodynamic theories so long as their energies are bounded below.
Abstract: We develop spatial-temporally second-order, energy stable numerical schemes for two classes of hydrodynamic phase field models of binary viscous fluid mixtures of different densities One is quasi-incompressible while the other is incompressible We introduce a novel energy quadratization technique to arrive at fully discrete linear schemes, where in each time step only a linear system needs to be solved These schemes are then shown to be unconditionally energy stable rigorously subject to periodic boundary conditions so that a large time step is plausible Both spatial and temporal mesh refinements are conducted to illustrate the second-order accuracy of the schemes The linearization technique developed in this paper is so general that it can be applied to any thermodynamically consistent hydrodynamic theories so long as their energies are bounded below Numerical examples on coarsening dynamics of two immiscible fluids and a heavy fluid drop settling in a lighter fluid matrix are presented to show the
78 citations
TL;DR: In this article, a fair performance comparison of the continuous finite element method, the symmetric interior penalty discontinuous Galerkin method, and the hybridized discontinuous GG method is presented.
Abstract: This study presents a fair performance comparison of the continuous finite element method, the symmetric interior penalty discontinuous Galerkin method, and the hybridized discontinuous Galerkin (H...
74 citations
TL;DR: A novel Hybrid High-Order method for the simulation of Darcy flows in fractured porous media that is fully robust with respect to the heterogeneity of the permeability coefficients, and it exhibits only a mild dependence on the square root of the local anisotropy of the bulk permeability.
Abstract: We develop a novel Hybrid High-Order method for the simulation of Darcy flows in fractured porous media. The discretization hinges on a mixed formulation in the bulk region and a primal formulation inside the fracture. Salient features of the method include a seamless treatment of nonconforming discretizations of the fracture, as well as the support of arbitrary approximation orders on fairly general meshes. For the version of the method corresponding to a polynomial degree $k\ge 0$, we prove convergence in $h^{k+1}$ of the discretization error measured in an energy-like norm. In the error estimate, we explicitly track the dependence of the constants on the problem data, showing that the method is fully robust with respect to the heterogeneity of the permeability coefficients, and it exhibits only a mild dependence on the square root of the local anisotropy of the bulk permeability. The numerical validation on a comprehensive set of test cases confirms the theoretical results.
74 citations
TL;DR: This paper considers the numerical approximation for a phase field model of the coupled two-phase free flow and two- phase porous media flow that consists of Cahn--Hilliard--Navier--Stokes equations in the free flow region and Cahn-Hilliard-Darcy-Navier-Stokes equation in the porous media region.
Abstract: In this paper, we consider the numerical approximation for a phase field model of the coupled two-phase free flow and two-phase porous media flow. This model consists of Cahn--Hilliard--Navier--Stokes equations in the free flow region and Cahn--Hilliard--Darcy equations in the porous media region that are coupled by seven interface conditions. The coupled system is decoupled based on the interface conditions and the solution values on the interface from the previous time step. A fully discretized scheme with finite elements for the spatial discretization is developed to solve the decoupled system. In order to deal with the difficulties arising from the interface conditions, the decoupled scheme needs to be constructed appropriately for the interface terms, and a modified discrete energy is introduced with an interface component. Furthermore, the scheme is linearized and energy stable. Hence, at each time step one need only solve a linear elliptic system for each of the two decoupled equations. Stability o...
TL;DR: In this article, the authors consider the discretization and subsequent model reduction of a system of partial differential-algebraic equations describing the propagation of pressure waves in a pipeline network.
Abstract: We consider the discretization and subsequent model reduction of a system of partial differential-algebraic equations describing the propagation of pressure waves in a pipeline network. Important properties like conservation of mass, dissipation of energy, passivity, existence of steady states, and exponential stability can be preserved by an appropriate semidiscretization in space via a mixed finite element method and also during the further dimension reduction by structure- preserving Galerkin projection, which is the main focus of this paper. Krylov subspace methods are employed for the construction of the reduced models, and we discuss certain modifications needed to satisfy some algebraic compatibility conditions which are required to ensure the well-posedness of the reduced models and the preservation of the key properties. Our arguments are based on a careful analysis of the underlying infinite dimensional problem and its Galerkin approximations. The proposed algorithms therefore have a direct inte...
TL;DR: In this paper, a second-order method for approximating the compressible Euler equations is introduced, which preserves all the known invariant domains of the Euler system: positivity of the density, posi...
Abstract: A new second-order method for approximating the compressible Euler equations is introduced. The method preserves all the known invariant domains of the Euler system: positivity of the density, posi...
TL;DR: A new approach for multi-material arbitrary Lagrangian--Eulerian (ALE) hydrodynamics simulations based on high-order finite elements posed onhigh-order curvilinear meshes, which depends critically on a functional perspective that enables subzonal physics and material modeling.
Abstract: We present a new approach for multi-material arbitrary Lagrangian--Eulerian (ALE) hydrodynamics simulations based on high-order finite elements posed on high-order curvilinear meshes. The method builds on and extends our previous work in the Lagrangian [V. A. Dobrev, T. V. Kolev, and R. N. Rieben, SIAM J. Sci. Comput., 34 (2012), pp. B606--B641] and remap [R. W. Anderson et al., Internat. J. Numer. Methods Fluids, 77 (2015), pp. 249--273] phases of ALE, and depends critically on a functional perspective that enables subzonal physics and material modeling [V. A. Dobrev et al., Internat. J. Numer. Methods Fluids, 82 (2016), pp. 689--706]. Curvilinear mesh relaxation is based on node movement, which is determined through the solution of an elliptic equation. The remap phase is posed in terms of advecting state variables between two meshes over a fictitious time interval. The resulting advection equation is solved by a discontinuous Galerkin (DG) formulation, combined with a customized Flux Corrected Transpor...
TL;DR: In this paper, a pseudospectral collocation approximation of the PDE dynamics and an iterative method for the nonlinear Hamilton-Jacobi-Bellman (HJB) equation associated to the feedback synthesis are proposed.
Abstract: A procedure for the numerical approximation of high-dimensional Hamilton--Jacobi--Bellman (HJB) equations associated to optimal feedback control problems for semilinear parabolic equations is proposed. Its main ingredients are a pseudospectral collocation approximation of the PDE dynamics and an iterative method for the nonlinear HJB equation associated to the feedback synthesis. The latter is known as the successive Galerkin approximation. It can also be interpreted as Newton iteration for the HJB equation. At every step, the associated linear generalized HJB equation is approximated via a separable polynomial approximation ansatz. Stabilizing feedback controls are obtained from solutions to the HJB equations for systems of dimension up to fourteen.
TL;DR: Numerical simulations in two and four dimensions for linear Landau damping and the two-stream instability highlight the favorable behavior of this numerical method and show that the proposed algorithm is able to drastically reduce the required computational effort.
Abstract: Many problems encountered in plasma physics require a description by kinetic equations, which are posed in an up to six-dimensional phase space. A direct discretization of this phase space, often c...
TL;DR: This paper proposes a block circulant preconditioner for the all-at-once evolutionary PDE system which has block Toeplitz structure, and rigorously establishes convergence bounds for MINRES which guarantee a number of iterations independent of the number of time-steps for theall- at-once system.
Abstract: Standard Krylov subspace solvers for self-adjoint problems have rigorous convergence bounds based solely on eigenvalues. However, for non-self-adjoint problems, eigenvalues do not determine behavior even for widely used iterative methods. In this paper, we discuss time-dependent PDE problems, which are always non-self-adjoint. We propose a block circulant preconditioner for the all-at-once evolutionary PDE system which has block Toeplitz structure. Through reordering of variables to obtain a symmetric system, we are able to rigorously establish convergence bounds for MINRES which guarantee a number of iterations independent of the number of time-steps for the all-at-once system. If the spatial differential operators are simultaneously diagonalizable, we are able to quickly apply the preconditioner through use of a sine transform; and for those that are not, we are able to use an algebraic multigrid process to provide a good approximation. Results are presented for solution to both the heat and convection ...
TL;DR: This work presents a systematic computational framework for generating positive quadrature rules in multiple dimensions on general geometries using penalty methods to address the geometric constraints, and subsequently solves a quadratic minimization problem via the Gauss-Newton method.
Abstract: We present a systematic computational framework for generating positive quadrature rules in multiple dimensions on general geometries. A direct moment-matching formulation that enforces exact integ...
TL;DR: TBLIS as mentioned in this paper implements tensor contraction using the flexible BLAS-like Instantiation Software (BLIS) framework, which allows transposition (reshaping) of the tensor to be fused with internal partitioning and packing operations, requiring no explicit transposition operations or additional workspace.
Abstract: Tensor computations---in particular tensor contraction (TC)---are important kernels in many scientific computing applications. Due to the fundamental similarity of TC to matrix multiplication and to the availability of optimized implementations such as the BLAS, tensor operations have traditionally been implemented in terms of BLAS operations, incurring both a performance and a storage overhead. Instead, we implement TC using the flexible BLAS-like Instantiation Software (BLIS) framework, which allows for transposition (reshaping) of the tensor to be fused with internal partitioning and packing operations, requiring no explicit transposition operations or additional workspace. This implementation, TBLIS, achieves performance approaching that of matrix multiplication, and in some cases considerably higher than that of traditional TC. Our implementation supports multithreading using an approach identical to that used for matrix multiplication in BLIS, with similar performance characteristics. The complexity...
TL;DR: This paper proposes comparing the methods of log-PCA and geodesic PCA in the Wasserstein space as introduced in [J. Bigot et al., Ann. Inst. Stat., 53 (2017)], and proposes a novel forward-backward algorithm that allows for a detailed comparison bet...
Abstract: This paper is concerned with the statistical analysis of datasets whose elements are random histograms. For the purpose of learning principal modes of variation from such data, we consider the issue of computing the principal component analysis (PCA) of histograms with respect to the 2-Wasserstein distance between probability measures. To this end, we propose comparing the methods of log-PCA and geodesic PCA in the Wasserstein space as introduced in [J. Bigot et al., Ann. Inst. Henri Poincare Probab. Stat., 53 (2017), pp. 1--26; V. Seguy and M. Cuturi, Principal geodesic analysis for probability measures under the optimal transport metric, in Advances in Neural Information Processing Systems 28, C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett, eds., Curran Associates, Inc., Red Hook, NY, 2015, pp. 3294--3302]. Geodesic PCA involves solving a nonconvex optimization problem. To solve it approximately, we propose a novel forward-backward algorithm. This allows us to give a detailed comparison bet...
TL;DR: The numerical approximations for a phase field model consisting of incompressible Navier--Stokes equations with a generalized Navier boundary condition, and the Cahn--Hilliard equation, are considered.
Abstract: We consider the numerical approximations for a phase field model consisting of incompressible Navier--Stokes equations with a generalized Navier boundary condition, and the Cahn--Hilliard equation ...
TL;DR: In this article, the authors consider a Stokes problem on a 2D surface embedded in a 3D domain and describe an equilibrium, area-preserving tangential flow of a viscous surface fluid.
Abstract: We consider a Stokes problem posed on a 2D surface embedded in a 3D domain. The equations describe an equilibrium, area-preserving tangential flow of a viscous surface fluid and serve as a model pr...
TL;DR: Improvements in the Firedrake finite element library are presented that allow for straightforward development of the building blocks of extensible, composable preconditioners that decouple the solver from the model formulation.
Abstract: The efficient solution of discretizations of coupled systems of partial differential equations (PDEs) is at the core of much of numerical simulation. Significant effort has been expended on scalable algorithms to precondition Krylov iterations for the linear systems that arise. With few exceptions, the reported numerical implementation of such solution strategies is specific to a particular model setup, and intimately ties the solver strategy to the discretization and PDE, especially when the preconditioner requires auxiliary operators. In this paper, we present recent improvements in the Firedrake finite element library that allow for straightforward development of the building blocks of extensible, composable preconditioners that decouple the solver from the model formulation. Our implementation extends the algebraic composability of linear solvers offered by the PETSc library by augmenting operators, and hence preconditioners, with the ability to provide any necessary auxiliary operators. Rather than s...
TL;DR: In this paper, the authors characterize a quadratic regularizer for graph transport with linear ground distance over a graph and theoretically analyze the behavior of quadratically regularized graph transport, characterizing how regularization affects the structure of flows in the regime of small but nonzero regularization.
Abstract: Optimal transportation provides a means of lifting distances between points on a geometric domain to distances between signals over the domain, expressed as probability distributions. On a graph, transportation problems can be used to express challenging tasks involving matching supply to demand with minimal shipment expense; in discrete language, these become minimum-cost network flow problems. Regularization typically is needed to ensure uniqueness for the linear ground distance case and to improve optimization convergence. In this paper, we characterize a quadratic regularizer for transport with linear ground distance over a graph. We theoretically analyze the behavior of quadratically regularized graph transport, characterizing how regularization affects the structure of flows in the regime of small but nonzero regularization. We further exploit elegant second-order structure in the dual of this problem to derive an easily implemented Newton-type optimization algorithm.
TL;DR: This paper designs a new type of high order finite volume weighted essentially nonoscillatory (WENO) schemes to solve hyperbolic conservation laws on triangular meshes and achieves any high order accuracy for the first time with the usage of only five unequal sized stencils.
Abstract: In this paper, we design a new type of high order finite volume weighted essentially nonoscillatory (WENO) schemes to solve hyperbolic conservation laws on triangular meshes. The main advantages of these schemes are their compactness and robustness and that they could maintain a good convergence property for some steady state problems. Compared with the classical finite volume WENO schemes [C. Hu and C.-W. Shu, J. Comput. Phys., 150 (1999), pp. 97--127], the optimal linear weights are independent of the topological structure of the triangular meshes and can be any positive numbers with the one requirement that their summation is one. This is the first time any high order accuracy with the usage of only five unequal sized stencils in a spatial reconstruction methodology on triangular meshes has been obtained. Extensive numerical results are provided to illustrate the good performance of such new finite volume WENO schemes.
TL;DR: A truncated Newton (TRN) method for time-domain full waveform inversion (FWI) in a visco-acoustic medium has been developed based on the 2nd-order adjoint state formulation.
Abstract: A truncated Newton (TRN) method for time-domain full waveform inversion (FWI) in a visco-acoustic medium has been developed based on the 2nd-order adjoint state formulation. Time-domain gradient es...
TL;DR: In this article, the authors present two different models to compute the pressure field and Darcy velocity in the system, one allows a normal flow out of a fracture at the intersections, while the second grants also a tangential flow along the intersections.
Abstract: Discrete fracture networks is a key ingredient in the simulation of physical processes which involve fluid flow in the underground, when the surrounding rock matrix is considered impervious. In this paper we present two different models to compute the pressure field and Darcy velocity in the system. The first allows a normal flow out of a fracture at the intersections, while the second grants also a tangential flow along the intersections. For the numerical discretization, we use the mixed virtual element method as it is known to handle grid elements of, almost, any arbitrary shape. The flexibility of the discretization allows us to loosen the requirements on grid construction, and thus significantly simplify the flow discretization compared to traditional discrete fracture network models. A coarsening algorithm, from the algebraic multigrid literature, is also considered to further speed up the computation. The performance of the method is validated by numerical experiments.
TL;DR: A randomized maximum a posteriori method for generating approximate samples of posteriors in high dimensional Bayesian inverse problems governed by large-scale forward problems and an approximate Metropolization approach is presented to reduce the bias in rMA.
Abstract: We present a randomized maximum a posteriori (rMAP) method for generating approximate samples of posteriors in high dimensional Bayesian inverse problems governed by large-scale forward problems. We derive the rMAP approach by (1) casting the problem of computing the MAP point as a stochastic optimization problem; (2) interchanging optimization and expectation; and (3) approximating the expectation with a Monte Carlo method. For a specific randomized data and prior mean, rMAP reduces to the randomized maximum likelihood (RML) approach. It can also be viewed as an iterative stochastic Newton method. An analysis of the convergence of the rMAP samples is carried out for both linear and nonlinear inverse problems. Each rMAP sample requires solution of a PDE-constrained optimization problem; to solve these problems, we employ a state-of-the-art trust region inexact Newton conjugate gradient method with sensitivity-based warm starts. An approximate Metropolization approach is presented to reduce the bias in rMA...
TL;DR: In this article, a new algorithm for automatic one-shot generation of scattered node sets on irregular two-dimensional and three-dimensional (3D) domains using Poisson disk sampling coupled to n...
Abstract: We present a new algorithm for the automatic one-shot generation of scattered node sets on irregular two dimensional (2D) and three dimensional (3D) domains using Poisson disk sampling coupled to n...
TL;DR: Second order, fully discrete, energy stable methods on spatially staggered grids for a hydrodynamic phase field model of binary viscous fluid mixtures in a confined geometry subject to both physical and periodic boundary conditions are presented.
Abstract: We present second order, fully discrete, energy stable methods on spatially staggered grids for a hydrodynamic phase field model of binary viscous fluid mixtures in a confined geometry subject to both physical and periodic boundary conditions. We apply the energy quadratization strategy to develop a linear-implicit scheme. We then extend it to a decoupled, linear scheme by introducing an intermediate velocity term so that the phase variable, velocity field, and pressure can be solved sequentially. The two new, fully discrete linear schemes are then shown to be unconditionally energy stable, and the linear systems resulting from the schemes are proved uniquely solvable. Rates of convergence of the two linear schemes in both space and time are verified numerically. The decoupled scheme tends to introduce excessive dissipation compared to the coupled one. The coupled scheme is then used to simulate fluid drops of one fluid in the matrix of another fluid as well as mixing dynamics of binary polymeric, viscous...