Do Current Lattice Boltzmann Methods for Diffusion and Advection-Diffusion Equations Respect Maximum Principle and the Non-Negative Constraint?

The Lattice Boltzmann Method (LBM) has established itself as a popular numerical method in computational fluid dynamics. Several advancements have been recently made in LBM, which include multiple-relaxation-time LBM to simulate anisotropic advection-diffusion processes. Because of the importance of LBM simulations for transport problems in subsurface and reactive flows, one needs to study the accuracy and structure preserving properties of numerical solutions under the LBM. The solutions to advective-diffusive systems are known to satisfy maximum principles, comparison principles, the non-negative constraint, and the decay property. In this paper, using several numerical experiments, it will be shown that current single- and multiple-relaxation-time lattice Boltzmann methods fail to preserve these mathematical properties for transient diffusion-type equations. We will also show that these violations may not be removed by simply refining the discretization parameters. More importantly, it will be shown that meeting stability conditions alone does not guarantee the preservation of the aforementioned mathematical principles and physical constraints in the discrete setting. A discussion on the source of these violations and possible approaches to avoid them is included. A condition to guarantee the non-negativity of concentration under LBM in the case of isotropic diffusion is also derived. The impact of this research is twofold. First, the study poses several outstanding research problems, which should guide researchers to develop LBM-based formulations for transport problems that respect important mathematical properties and physical constraints in the discrete setting. This paper can also serve as a good source of benchmark problems for such future research endeavors. Second, this study cautions the practitioners of the LBM for transport problems with the associated numerical deficiencies of the LBM, and provides guidelines for performing predictive simulations of advective-diffusive processes using the LBM.     

Mechanics-Based Solution Verification for Porous Media Models

This paper presents a new approach to verify the accuracy of computational simulations. We develop mathematical theorems which can serve as robust a posteriori error estimation techniques to identify numerical pollution, check the performance of adaptive meshes, and verify numerical solutions. We demonstrate performance of this methodology on problems from flow thorough porous media. However, one can extend it to other models. We construct mathematical properties such that the solutions to Darcy and Darcy-Brinkman equations satisfy them. The mathematical properties include the total minimum mechanical power, minimum dissipation theorem, reciprocal relation, and maximum principle for the vorticity. All the developed theorems have firm mechanical bases and are independent of numerical methods. So, these can be utilized for solution verification of finite element, finite volume, finite difference, lattice Boltzmann methods and so forth. In particular, we show that, for a given set of boundary conditions, Darcy velocity has the minimum total mechanical power of all the kinematically admissible vector fields. We also show that a similar result holds for Darcy-Brinkman velocity. We then show for a conservative body force, the Darcy and Darcy-Brinkman velocities have the minimum total dissipation among their respective kinematically admissible vector fields. Using numerical examples, we show that the minimum dissipation and total mechanical power theorems can be utilized to identify pollution errors in numerical solutions. The solutions to Darcy and Darcy-Brinkman equations are shown to satisfy a reciprocal relation, which has the potential to identify errors in the numerical implementation of boundary conditions. It is also shown that the vorticity under both steady and transient Darcy-Brinkman equations satisfy maximum principles if the body force is conservative and the permeability is homogeneous and isotropic. A discussion on the nature of vorticity under steady and transient Darcy equations is also presented. Using several numerical examples, we will demonstrate the predictive capabilities of the proposed a posteriori techniques in assessing the accuracy of numerical solutions for a general class of problems, which could involve complex domains and general computational grids.     

Constraint Preserving Schemes Using Potential-Based Fluxes I Multidimensional Transport Equations

We consider constraint preserving multidimensional evolution equations. A prototypical example is provided by the magnetic induction equation of plasma physics. The constraint of interest is the divergence of the magnetic field. We design finite volume schemes which approximate these equations in a stable manner and preserve a discrete version of the constraint. The schemes are based on reformulating standard edge centered finite volume fluxes in terms of vertex centered potentials. The potential-based approach provides a general framework for faithful discretizations of constraint transport and we apply it to both divergence preserving as well as curl preserving equations. We present benchmark numerical tests which confirm that our potential-based schemes achieve high resolution, while being constraint preserving.     

Truncation Errors,Exact and Heuristic Stability Analysis of Two-Relaxation-Times Lattice Boltzmann Schemes for Anisotropic Advection-Diffusion Equation

This paper establishes relations between the stability and the high-order truncated corrections for modeling of the mass conservation equation with the two-relaxation-times (TRT) collision operator. First we propose a simple method to derive the truncation errors from the exact, central-difference type, recurrence equations of the TRT scheme. They also supply its equivalent three-time-level discretization form. Two different relationships of the two relaxation rates nullify the third (advection) and fourth (pure diffusion) truncation errors, for any linear equilibrium and any velocity set. However, the two relaxation times alone cannot remove the leading-order advection-diffusion error, because of the intrinsic fourth-order numerical diffusion. The truncation analysis is carefully verified for the evolution of concentration waves with the anisotropic diffusion tensors. The anisotropic equilibrium functions are presented in a simple but general form, suitable for the minimal velocity sets and the d2Q9, d3Q13, d3Q15 and d3Q19 velocity sets. All anisotropic schemes are complemented by their exact necessary von Neumann stability conditions and equivalent finite-difference stencils. The sufficient stability conditions are proposed for the most stable (OTRT) family, which enables modeling at any Peclet numbers with the same velocity amplitude. The heuristic stability analysis of the fourth-order truncated corrections extends the optimal stability to larger relationships of the two relaxation rates, in agreement with the exact (one-dimensional) and numerical (multi-dimensional) stability analysis. A special attention is put on the choice of the equilibrium weights. By combining accuracy and stability predictions, several strategies for selecting the relaxation and free-tunable equilibrium parameters are suggested and applied to the evolution of the Gaussian hill.     

Discrete Maximum Principle for the Weak Galerkin Method for Anisotropic Diffusion Problems

A weak Galerkin discretization of the boundary value problem of a general anisotropic diffusion problem is studied for preservation of the maximum principle. It is shown that the direct application of the M-matrix theory to the stiffness matrix of the weak Galerkin discretization leads to a strong mesh condition requiring all of the mesh dihedral angles to be strictly acute (a constant-order away from 90 degrees). To avoid this difficulty, a reduced system is considered and shown to satisfy the discrete maximum principle under weaker mesh conditions. The discrete maximum principle is then established for the full weak Galerkin approximation using the relations between the degrees of freedom located on elements and edges. Sufficient mesh conditions for both piecewise constant and general anisotropic diffusion matrices are obtained. These conditions provide a guideline for practical mesh generation for preservation of the maximum principle. Numerical examples are presented.     

A Modeling Framework for Coupling Plasticity with Species Diffusion

This paper presents a modeling framework—mathematical model and computational framework—to study the response of a plastic material due to the presence and transport of a chemical species in the host material. Such a modeling framework is important to a wide variety of problems ranging from Li-ion batteries, moisture diffusion in cementitious materials, hydrogen diffusion in metals, to consolidation of soils under severe loading-unloading regimes. The mathematical model incorporates experimental observations reported in the literature on how (elastic and plastic) material properties change because of the presence and transport of a chemical species. Also, the model accounts for one-way (transport affects the deformation but not vice versa) and two-way couplings between deformation and transport subproblems. The resulting coupled equations are not amenable to analytical solutions; so, we present a robust computational framework for obtaining numerical solutions. Given that popular numerical formulations do not produce nonnegative solutions, the computational framework uses an optimized-based nonnegative formulation that respects physical constraints (e.g., nonnegative concentrations). For completeness, we also show the effect and propagation of the negative concentrations, often produced by contemporary transport solvers, into the overall predictions of deformation and concentration fields. Notably, anisotropy of the diffusion process exacerbates these unphysical violations. Using representative numerical examples, we discuss how the concentration field affects plastic deformations of a degrading solid. Based on these numerical examples, we also discuss how plastic zones spread because of material degradation. To illustrate how the proposed computational framework performs, we report various performance metrics such as optimization iterations and time-to-solution.     

A Weighted Runge-Kutta Discontinuous Galerkin Method for 3D Acoustic and Elastic Wave-Field Modeling

Numerically solving 3D seismic wave equations is a key requirement for forward modeling and inversion. Here, we propose a weighted Runge-Kutta discontinuous Galerkin (WRKDG) method for 3D acoustic and elastic wave-field modeling. For this method, the second-order seismic wave equations in 3D heterogeneous anisotropic media are transformed into a first-order hyperbolic system, and then we use a discontinuous Galerkin (DG) solver based on numerical-flux formulations for spatial discretization. The time discretization is based on an implicit diagonal Runge-Kutta (RK) method and an explicit iterative technique, which avoids solving a large-scale system of linear equations. In the iterative process, we introduce a weighting factor. We investigate the numerical stability criteria of the 3D method in detail for linear and quadratic spatial basis functions. We also present a 3D analysis of numerical dispersion for the full discrete approximation of acoustic equation, which demonstrates that the WRKDG method can efficiently suppress numerical dispersion on coarse grids. Numerical results for several different 3D models including homogeneous and heterogeneous media with isotropic and anisotropic cases show that the 3D WRKDG method can effectively suppress numerical dispersion and provide accurate wave-field information on coarse mesh.     

An Energy Regularization Method for the Backward Diffusion Problem and Its Applications to Image Deblurring

For the backward diffusion equation, a stable discrete energy regularization algorithm is proposed. Existence and uniqueness of the numerical solution are given. Moreover, the error between the solution of the given backward diffusion equation and the numerical solution via the regularization method can be estimated. Some numerical experiments illustrate the efficiency of the method, and its application in image deblurring.     

A Nonlinear Finite Volume Scheme Preserving Maximum Principle for Diffusion Equations

In this paper we propose a new nonlinear cell-centered finite volume scheme on general polygonal meshes for two dimensional anisotropic diffusion problems, which preserves discrete maximum principle (DMP). The scheme is based on the so-called diamond scheme with a nonlinear treatment on its tangential flux to obtain a local maximum principle (LMP) structure. It is well-known that existing DMP preserving diffusion schemes suffer from the fact that auxiliary unknowns should be presented as a convex combination of primary unknowns. In this paper, to get rid of this constraint a nonlinearization strategy is introduced and it requires only a second-order accurate approximation for auxiliary unknowns. Numerical results show that this scheme has second-order accuracy, preserves maximum and minimum for solutions and is conservative.     

Enforcing the Discrete Maximum Principle for Linear Finite Element Solutions of Second-Order Elliptic Problems

The maximum principle is a basic qualitative property of the solution of second-order elliptic boundary value problems. The preservation of the qualitative characteristics, such as the maximum principle, in discrete model is one of the key requirements. It is well known that standard linear finite element solution does not satisfy maximum principle on general triangular meshes in 2D. In this paper we consider how to enforce discrete maximum principle for linear finite element solutions for the linear second-order self-adjoint elliptic equation. First approach is based on repair technique, which is a posteriori correction of the discrete solution. Second method is based on constrained optimization. Numerical tests that include anisotropic cases demonstrate how our method works for problems for which the standard finite element methods produce numerical solutions that violate the discrete maximum principle.     

A Nearly Analytical Discrete Method for Wave-Field Simulations in 2D Porous Media

The nearly analytic discrete method (NADM) is a perturbation method originally proposed by Yang et al. (2003) [26] for acoustic and elastic waves in elastic media. This method is based on a truncated Taylor series expansion and interpolation approximations and it can suppress effectively numerical dispersions caused by the discretizating the wave equations when too-coarse grids are used. In the present work, we apply the NADM to simulating acoustic and elastic wave propagations in 2D porous media. Our method enables wave propagation to be simulated in 2D porous isotropic and anisotropic media. Numerical experiments show that the error of the NADM for the porous case is less than those of the conventional finite-difference method (FDM) and the so-called Lax-Wendroff correction (LWC) schemes. The three-component seismic wave fields in the 2D porous isotropic medium are simulated and compared with those obtained by using the LWC method and exact solutions. Several characteristics of wave propagating in porous anisotropic media, computed by the NADM, are also reported in this study. Promising numerical results illustrate that the NADM provides a useful tool for large-scale porous problems and it can suppress effectively numerical dispersions.     

Extrapolation Cascadic Multigrid Method for Cell-Centered FV Discretization of Diffusion Equations with Strongly Discontinuous and Anisotropic Coefficients

Extrapolation cascadic multigrid (EXCMG) method with conjugate gradient smoother is very efficient for solving the elliptic boundary value problems with linear finite element discretization. However, it is not trivial to generalize the vertex-centred EXCMG method to cell-centered finite volume (FV) methods for diffusion equations with strongly discontinuous and anisotropic coefficients, since a non-nested hierarchy of grid nodes are used in the cell-centered discretization. For cell-centered FV schemes, the vertex values (auxiliary unknowns) need to be approximated by cell-centered ones (primary unknowns). One of the novelties is to propose a new gradient transfer (GT) method of interpolating vertex unknowns with cell-centered ones, which is easy to implement and applicable to general diffusion tensors. The main novelty of this paper is to design a multigrid prolongation operator based on the GT method and splitting extrapolation method, and then propose a cell-centered EXCMG method with BiCGStab smoother for solving the large linear system resulting from linear FV discretization of diffusion equations with strongly discontinuous and anisotropic coefficients. Numerical experiments are presented to demonstrate the high efficiency of the proposed method.     

Efficient Energy Stable Schemes with Spectral Discretization in Space for Anisotropic Cahn-Hilliard Systems

We develop in this paper efficient and robust numerical methods for solving anisotropic Cahn-Hilliard systems. We construct energy stable schemes for the time discretization of the highly nonlinear anisotropic Cahn-Hilliard systems by using a stabilization technique. At each time step, these schemes lead to a sequence of linear coupled elliptic equations with constant coefficients that can be efficiently solved by using a spectral-Galerkin method. We present numerical results that are consistent with earlier work on this topic, and also carry out various simulations, such as the linear bi-Laplacian regularization and the nonlinear Willmore regularization, to demonstrate the efficiency and robustness of the new schemes.     

A High Order Sharp-Interface Method with Local Time Stepping for Compressible Multiphase Flows

In this paper, a new sharp-interface approach to simulate compressible multiphase flows is proposed. The new scheme consists of a high order WENO finite volume scheme for solving the Euler equations coupled with a high order path-conservative discontinuous Galerkin finite element scheme to evolve an indicator function that tracks the material interface. At the interface our method applies ghost cells to compute the numerical flux, as the ghost fluid method. However, unlike the original ghost fluid scheme of Fedkiw et al. [15], the state of the ghost fluid is derived from an approximate-state Riemann solver, similar to the approach proposed in [25], but based on a much simpler formulation. Our formulation leads only to one single scalar nonlinear algebraic equation that has to be solved at the interface, instead of the system used in [25]. Away from the interface, we use the new general Osher-type flux recently proposed by Dumbser and Toro [13], which is a simple but complete Riemann solver, applicable to general hyperbolic conservation laws. The time integration is performed using a fully-discrete one-step scheme, based on the approaches recently proposed in [5, 7]. This allows us to evolve the system also with time-accurate local time stepping. Due to the sub-cell resolution and the subsequent more restrictive time-step constraint of the DG scheme, a local evolution for the indicator function is applied, which is matched with the finite volume scheme for the solution of the Euler equations that runs with a larger time step. The use of a locally optimal time step avoids the introduction of excessive numerical diffusion in the finite volume scheme. Two different fluids have been used, namely an ideal gas and a weakly compressible fluid modeled by the Tait equation. Several tests have been computed to assess the accuracy and the performance of the new high order scheme. A verification of our algorithm has been carefully carried out using exact solutions as well as a comparison with other numerical reference solutions. The material interface is resolved sharply and accurately without spurious oscillations in the pressure field.     

A Numerical Method and Its Error Estimates for the Decoupled Forward-Backward Stochastic Differential Equations

In this paper, a new numerical method for solving the decoupled forward-backward stochastic differential equations (FBSDEs) is proposed based on some specially derived reference equations. We rigorously analyze errors of the proposed method under general situations. Then we present error estimates for each of the specific cases when some classical numerical schemes for solving the forward SDE are taken in the method; in particular, we prove that the proposed method is second-order accurate if used together with the order-2.0 weak Taylor scheme for the SDE. Some examples are also given to numerically demonstrate the accuracy of the proposed method and verify the theoretical results.     

Quantum Implementation of Numerical Methods for Convection-Diffusion Equations: Toward Computational Fluid Dynamics

We present quantum numerical methods for the typical initial boundary value problems (IBVPs) of convection-diffusion equations in fluid dynamics. The IBVP is discretized into a series of linear systems via finite difference methods and explicit time marching schemes. To solve these discrete systems in quantum computers, we design a series of quantum circuits, including four stages of encoding, amplification, adding source terms, and incorporating boundary conditions. In the encoding stage, the initial condition is encoded in the amplitudes of quantum registers as a state vector to take advantage of quantum algorithms in space complexity. In the following three stages, the discrete differential operators in classical computing are converted into unitary evolutions to satisfy the postulate in quantum systems. The related arithmetic calculations in quantum amplitudes are also realized to sum up the increments from these stages. The proposed quantum algorithm is implemented within the open-source quantum computing framework Qiskit [2]. By simulating one-dimensional transient problems, including the Helmholtz equation, the Burgers’ equation, and Navier-Stokes equations, we demonstrate the capability of quantum computers in fluid dynamics.     

A Positivity-Preserving Second-Order BDF Scheme for the Cahn-Hilliard Equation with Variable Interfacial Parameters

We present and analyze a new second-order finite difference scheme for the Macromolecular Microsphere Composite hydrogel, Time-Dependent Ginzburg-Landau (MMC-TDGL) equation, a Cahn-Hilliard equation with Flory-Huggins-deGennes energy potential. This numerical scheme with unconditional energy stability is based on the Backward Differentiation Formula (BDF) method in time derivation combining with Douglas-Dupont regularization term. In addition, we present a pointwise bound of the numerical solution for the proposed scheme in the theoretical level. For the convergent analysis, we treat three nonlinear logarithmic terms as a whole and deal with all logarithmic terms directly by using the property that the nonlinear error inner product is always non-negative. Moreover, we present the detailed convergent analysis in $ℓ^∞$(0,$T$;$H_h^{-1}$)∩$ℓ^2$(0,$T$;$H_h^1$) norm. At last, we use the local Newton approximation and multigrid method to solve the nonlinear numerical scheme, and various numerical results are presented, including the numerical convergence test, positivity-preserving property test, spinodal decomposition, energy dissipation and mass conservation properties.     

A Riccati-equation-based algorithm for continuous-time optimal control problems

In this paper we consider continuous-time unconstrained optimal control problems. We propose a computational method which is essentially based on the closed-loop solutions of the linear quadratic optimal control problems. In the proposed algorithm, Riccati differential equations play an important role. We prove that accumulation points generated by the present algorithm, if they exist, satisfy the weak necessary conditions for optimality, under some assumptions including Kalman's sufficient conditions for the bounded Riccati solutions. In addition, we also propose the simple but effective technique to guarantee the boundedness of the solutions of Riccati equations. Lastly, we illustrate the usefulness of the present algorithm through simulation experiences. Copyright © 1998 John Wiley & Sons, Ltd.     

On Universal Osher-Type Schemes for General Nonlinear Hyperbolic Conservation Laws

This paper is concerned with a new version of the Osher-Solomon Riemann solver and is based on a numerical integration of the path-dependent dissipation matrix. The resulting scheme is much simpler than the original one and is applicable to general hyperbolic conservation laws, while retaining the attractive features of the original solver: the method is entropy-satisfying, differentiable and complete in the sense that it attributes a different numerical viscosity to each characteristic field, in particular to the intermediate ones, since the full eigenstructure of the underlying hyperbolic system is used. To illustrate the potential of the proposed scheme we show applications to the following hyperbolic conservation laws: Euler equations of compressible gasdynamics with ideal gas and real gas equation of state, classical and relativistic MHD equations as well as the equations of nonlinear elasticity. To the knowledge of the authors, apart from the Euler equations with ideal gas, an Osher-type scheme has never been devised before for any of these complicated PDE systems. Since our new general Riemann solver can be directly used as a building block of high order finite volume and discontinuous Galerkin schemes we also show the extension to higher order of accuracy and multiple space dimensions in the new framework of P_NP_M schemes on unstructured meshes recently proposed in [9].     

A Well-Posed and Discretely Stable Perfectly Matched Layer for Elastic Wave Equations in Second Order Formulation

We present a well-posed and discretely stable perfectly matched layer for the anisotropic (and isotropic) elastic wave equations without first re-writing the governing equations as a first order system. The new model is derived by the complex coordinate stretching technique. Using standard perturbation methods we show that complex frequency shift together with a chosen real scaling factor ensures the decay of eigen-modes for all relevant frequencies. To buttress the stability properties and the robustness of the proposed model, numerical experiments are presented for anisotropic elastic wave equations. The model is approximated with a stable node-centered finite difference scheme that is second order accurate both in time and space.