Positivity Preserving Research Articles

Efficient numerical solution of the Fokker–Planck equation using physics-conforming finite element methods

Abstract We consider the Fokker–Planck equation (FPE) for the orientation probability density of fiber suspensions. Using the continuous Galerkin method, we express the numerical solution in terms of Lagrange basis functions that are associated with N nodes of a computational mesh for a domain in the 3D physical space and M nodes of a mesh for the surface of a unit sphere representing the configuration space. The NM time-dependent unknowns of our finite element approximations are probabilities corresponding to discrete space locations and orientation angles. The framework of alternating-direction methods enables us to update the numerical solution in parallel by solving N evolution equations on the sphere and M three-dimensional advection equations in each (pseudo-)time step. To ensure positivity preservation as well as the normalization property of the probability density, we perform algebraic flux correction for each equation and synchronize the correction factors corresponding to different orientation angles. The velocity field for the spatial advection step is obtained using a Schur complement method to solve a generalized system of the incompressible Navier–Stokes equations (NSE). Fiber-induced subgrid-scale effects are taken into account using an effective stress tensor that depends on the second- and fourth-order moments of the orientation density function. Numerical studies are performed for individual subproblems and for the coupled FPE-NSE system.

On high order positivity-preserving well-balanced finite volume methods for the Euler equations with gravitation

This paper studies three high-order structure-preserving finite volume weighted essentially non-oscillatory (WENO) methods, which are not only well balanced (WB) for a general known hydrostatic equilibrium state but also preserve the positivity of density and pressure, for the compressible Euler equations under gravitational fields. These methods are built on a simple local scaling positivity-preserving (PP) limiter and a modified WENO-ZQ reconstruction exactly preserving the cell average value and scaling invariance. The WB properties of these three methods are achieved based on suitable numerical fluxes and approximation to the gravitational source terms. Based on some convex decomposition techniques as well as several critical properties of the admissible states and numerical flux, we carry out rigorous positivity-preserving analyses for these three WB schemes. We rigorously prove that the three WB methods, coupled with the PP limiter and a strong-stability-preserving time discretization, are always PP under suitable Courant-Friedrichs-Lewy conditions. Extensive numerical examples are provided to confirm WB and PP properties of three methods.

A Comparative Study of Different Sets of Variables in a Discontinuous Galerkin Method with Entropy Balance Enforcement

This paper investigates the effectiveness of different sets of variables in solving the compressible Euler equations using a modal Discontinuous Galerkin framework. Alongside the commonly used conservative and primitive variables, the entropy and logarithmic sets are considered to enforce entropy conservation/stability and positivity preservation of the thermodynamic state, respectively. An explicit correction to enforce entropy conservation/stability at the discrete level is also considered, with a significant increase in robustness for some of the solution strategies. Several two-dimensional inviscid test cases are computed to compare the performance of the different sets of variables, adding a directional shock-capturing term to the discretised equations when necessary. The entropy and logarithmic sets proved to be the most robust, completing simulations of an astrophysical jet at Mach number 2000 up to polynomial degree seven.

A positivity preserving high‐order finite difference method for compressible two‐fluid flows

AbstractAny robust computational scheme for compressible flows must retain the hyperbolicity property or the real‐valued sound speed. Failure to maintain hyperbolicity, or the positivity of the square of the speed of sound, causes nonphysical distortions and the blow‐up of numerical simulations. Strong shock waves and interfacial discontinuities are ubiquitous features of the two‐fluid compressible dynamics that can potentially induce positivity‐related failure in a simulation. This article presents a positivity‐preserving algorithm for a high‐order, primitive variable‐based, weighted essentially non‐oscillatory finite difference scheme. The positivity preservation relies on a flux limiting technique that locally adapts high‐order fluxes towards the first order to retain the physical bounds of the solution without loss of high‐order convergence. This positivity preserving scheme has been devised and implemented up to eleventh order in one and two dimensions for a two‐fluid compressible model that consists of a single mass, momentum, and energy equations, as well as an advection of material parameters for capturing the interfaces. Several one‐ and two‐dimensional challenging test problems verify the performance. The scheme effectively retains high order accuracy while allowing for the simulation of several challenging problems that otherwise could not be successfully solved using the base scheme, without any penalty on the CFL condition requirement, and without any significant impact on the CPU times. The scheme represents the first high‐order (up to 11‐order) hyperbolicity‐preserving scheme for the considered two‐fluid compressible flows in the fully Eulerian formulation. The inherent efficiency of finite differences and the new robust positivity preserving quality enable modeling other challenging problems of two‐fluid and two‐phase problems.

A Well-Balanced Active Flux Method for the Shallow Water Equations with Wetting and Drying

Active Flux is a third order accurate numerical method which evolves cell averages and point values at cell interfaces independently. It naturally uses a continuous reconstruction, but is stable when applied to hyperbolic problems. In this work, the Active Flux method is extended for the first time to a nonlinear hyperbolic system of balance laws, namely, to the shallow water equations with bottom topography. We demonstrate how to achieve an Active Flux method that is well-balanced, positivity preserving, and allows for dry states in one spatial dimension. Because of the continuous reconstruction all these properties are achieved using new approaches. To maintain third order accuracy, we also propose a novel high-order approximate evolution operator for the update of the point values. A variety of test problems demonstrates the good performance of the method even in presence of shocks.

Some standard and nonstandard finite difference schemes for a reaction–diffusion–chemotaxis model

Abstract Two standard and two nonstandard finite difference schemes are constructed to solve a basic reaction–diffusion–chemotaxis model, for which no exact solution is known. The continuous model involves a system of nonlinear coupled partial differential equations subject to some specified initial and boundary conditions. It is not possible to obtain theoretically the stability region of the two standard finite difference schemes. Through running some numerical experiments, we deduce heuristically that these classical methods give reasonable solutions when the temporal step size k k is chosen such that k ≤ 0.25 k\le 0.25 with the spatial step size h h fixed at h = 1.0 h=1.0 (first novelty of this work). We observe that the standard finite difference schemes are not always positivity preserving, and this is why we consider nonstandard finite difference schemes. Two nonstandard methods abbreviated as NSFD1 and NSFD2 from Chapwanya et al. are considered. NSFD1 was not used by Chapwanya et al. to generate results for the basic reaction–diffusion–chemotaxis model. We find that NSFD1 preserves positivity of the continuous model if some criteria are satisfied, namely, ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 γ ≤ 1 2 σ + β \frac{\phi \left(k)}{{\left[\psi \left(h)]}^{2}}=\frac{1}{2\gamma }\le \frac{1}{2\sigma +\beta } and β ≤ σ \beta \le \sigma , and this is the second novelty of this work. Chapwanya et al. modified NSFD1 to obtain NSFD2, which is positivity preserving if R = ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 γ R=\frac{\phi \left(k)}{{\left[\psi \left(h)]}^{2}}=\frac{1}{2\gamma } and 2 σ R ≤ 1 2\sigma R\le 1 , that is σ ≤ γ \sigma \le \gamma , and they presented some results. For the third highlight of this work, we show that NSFD2 is not always consistent and prove that consistency can be achieved if β → 0 \beta \to 0 and k h 2 → 0 \frac{k}{{h}^{2}}\to 0 . Fourthly, we show numerically that the rate of convergence in time of the four methods for case 2 is approximately one.

Provably Positive Central Discontinuous Galerkin Schemes via Geometric Quasilinearization for Ideal MHD Equations

In the numerical simulation of ideal magnetohydrodynamics (MHD), keeping the pressure and density always positive is essential for both physical considerations and numerical stability. This is however a challenging task, due to the underlying relation between such a positivity-preserving (PP) property and the magnetic divergence-free (DF) constraint as well as the strong nonlinearity of the MHD equations. In this paper, we present the first rigorous PP analysis of the central discontinuous Galerkin (CDG) methods and construct arbitrarily high-order provably PP CDG schemes for ideal MHD. By the recently developed geometric quasilinearization (GQL) approach, our analysis reveals that the PP property of standard CDG methods is closely related to a discrete magnetic DF condition, whose form was unknown prior to our analysis and differs from that for the noncentral discontinuous Galerkin (DG) and finite volume methods in [K. Wu, SIAM J. Numer. Anal., 56 (2018), pp. 2124–2147]. The discovery of this relation lays the foundation for the design of our PP CDG schemes. In the one-dimensional (1D) case, the discrete DF condition is naturally satisfied, and we rigorously prove that the standard CDG method is PP under a condition that can be enforced easily with an existing PP limiter. However, in the multidimensional cases, the corresponding discrete DF condition is highly nontrivial, yet critical, and we analytically prove that the standard CDG method, even with the PP limiter, is not PP in general, as it generally fails to meet the discrete DF condition. We address this issue by carefully analyzing the structure of the discrete divergence terms and then constructing new locally DF CDG schemes for Godunov’s modified MHD equations with an additional source term. The key point is to find out the suitable discretization of the source term such that it exactly cancels out all the terms in the discovered discrete DF condition. Based on the GQL approach, we prove in theory the PP property of the new multidimensional CDG schemes under a CFL condition. The robustness and accuracy of the proposed PP CDG schemes are further validated by several demanding 1D, two-dimensional, and three-dimensional numerical MHD examples, including the high-speed jets and blast problems with very low plasma beta.

A diffuse-domain-based numerical method for a chemotaxis-fluid model

In this paper, we consider a coupled chemotaxis-fluid system that models self-organized collective behavior of oxytactic bacteria in a sessile drop. This model describes the biological chemotaxis phenomenon in the fluid environment and couples a convective chemotaxis system for the oxygen-consuming and oxytactic bacteria with the incompressible Navier–Stokes equations subject to a gravitational force, which is proportional to the relative surplus of the cell density compared to the water density. We develop a new positivity preserving and high-resolution method for the studied chemotaxis-fluid system. Our method is based on the diffuse-domain approach, which we use to derive a new chemotaxis-fluid diffuse-domain (cf-DD) model for simulating bioconvection in complex geometries. The drop domain is imbedded into a larger rectangular domain, and the original boundary is replaced by a diffuse interface with finite thickness. The original chemotaxis-fluid system is reformulated on the larger domain with additional source terms that approximate the boundary conditions on the physical interface. We show that the cf-DD model converges to the chemotaxis-fluid model asymptotically as the width of the diffuse interface shrinks to zero. We numerically solve the resulting cf-DD system by a second-order hybrid finite-volume finite-difference method and demonstrate the performance of the proposed approach on a number of numerical experiments that showcase several interesting chemotactic phenomena in sessile drops of different shapes, where the bacterial patterns depend on the droplet geometries.

Bounds preserving temporal integration methods for hyperbolic conservation laws

In this work, we present a modification of explicit Runge–Kutta temporal integration schemes that guarantees the preservation of any locally-defined quasiconvex set of bounds for the solution. These schemes operate on the basis of a bijective mapping between an admissible set of solutions and the real domain to strictly enforce bounds. Within this framework, we show that it is possible to recover a wide range of methods independently of the spatial discretization, including positivity preserving, discrete maximum principle satisfying, entropy dissipative, and invariant domain preserving schemes. Furthermore, these schemes are proven to recover the order of accuracy of the underlying Runge–Kutta method upon which they are built. The additional computational cost is the evaluation of two nonlinear mappings which generally have closed-form solutions. We show the utility of this approach in numerical experiments using a pseudospectral spatial discretization without any explicit shock capturing schemes for nonlinear hyperbolic problems with discontinuities.

“A posteriori” limited high order and robust schemes for transient simulations of fluid flows in gas dynamics

In this paper, we propose a novel approximation strategy for time-dependent hyperbolic systems of conservation laws for the Euler system of gas dynamics that aims to represent the dynamics of strong interacting discontinuities. The goal of our method is to allow an approximation with a high-order of accuracy in smooth regions of the flow, while ensuring robustness and a non-oscillatory behaviour in the regions of steep gradients, in particular across shocks.Following the Multidimensional Optimal Order Detection (MOOD) ([15,17]) approach, a candidate solution is computed at a next time level via a high-order accurate explicit scheme ([3,5]). A so-called detector determines if the candidate solution reveals any spurious oscillations or numerical issue and, if so, only the troubled cells are locally recomputed via a more dissipative scheme. This allows to design a family of “a posteriori” limited, robust and positivity preserving, as well as high accurate, non-oscillatory and effective scheme. Among the detecting criteria of the novel MOOD strategy, two different approaches from literature, based on the work of [15,17] and on [35], are investigated. Numerical examples in 1D and 2D, on structured and unstructured meshes, are proposed to assess the effective order of accuracy for smooth flows, the non-oscillatory behaviour on shocked flows, the robustness and positivity preservation on more extreme flows.

Fifth-order well-balanced positivity-preserving finite difference AWENO scheme with hydrostatic reconstruction for hyperbolic chemotaxis models

We propose a fifth-order well-balanced positivity-preserving finite difference alternative WENO (AWENO) scheme with the affine-invariant WENO interpolation based on the Z-type nonlinear weights for the hyperbolic chemotaxis models. By using the techniques of source term reformulation, hydrostatic reconstruction of the interpolated conservative variables and modification of the numerical fluxes, the finite difference discretization is fifth-order and well-balanced. Moreover, the first-order interpolation with the Lax-Friedrichs (LF) flux and a reduced time step for the proposed discretized scheme has been shown to satisfy the density's positivity-preserving (PP) property. Thus, a simple positivity-preserving (PP) limiter conjugating the fifth-order hydrostatic reconstructed flux with the first-order positivity-preserving LF flux is introduced for extreme problems. Meanwhile, this improved approach strictly guarantees well-balanced property at the discrete level. Finally, one-, two-, and three-dimensional numerical examples are given to demonstrate the performance of the proposed AWENO scheme for this class of chemotaxis problems.

A nonstandard finite difference scheme for a time-fractional model of Zika virus transmission.

In this work, we investigate the transmission dynamics of the Zika virus, considering both a compartmental model involving humans and mosquitoes and an extended model that introduces a non-human primate (monkey) as a second reservoir host. The novelty of our approach lies in the later generalization of the model using a fractional time derivative. The significance of this study is underscored by its contribution to understanding the complex dynamics of Zika virus transmission. Unlike previous studies, we incorporate a non-human primate reservoir host into the model, providing a more comprehensive representation of the disease spread. Our results reveal the importance of utilizing a nonstandard finite difference (NSFD) scheme to simulate the disease's dynamics accurately. This NSFD scheme ensures the positivity of the solution and captures the correct asymptotic behavior, addressing a crucial limitation of standard solvers like the Runge-Kutta Fehlberg method (ode45). The numerical simulations vividly demonstrate the advantages of our approach, particularly in terms of positivity preservation, offering a more reliable depiction of Zika virus transmission dynamics. From these findings, we draw the conclusion that considering a non-human primate reservoir host and employing an NSFD scheme significantly enhances the accuracy and reliability of modeling Zika virus transmission. Researchers and policymakers can use these insights to develop more effective strategies for disease control and prevention.

A Novel Sediment Transport Model (STM) Accounting Phase Lag Effect. A Resonance Condition

The classical Exner model coupled with a bed-load sediment flux formula is widely used to describe the morphodynamics of coastal environments. However, the main drawbacks of this model are (i) Lack of robustness, (ii) Lack of differentiation between sediment and fluid velocities, and (iii) Generation of instabilities when the interactions between sediment and fluid flow become more important. Moreover, Exner's model does not allow us to know with which characteristic velocity the bottom is moving. This set of drawbacks weakens the effectiveness of most sediment transport models proposed in the literature, particularly the Exner model. In this work, we reformulate the bed-load equation and we propose a new averaged sediment transport model for application in coastal or estuarine environments. The proposed model incorporates phase shift effects into the bed-load equation. The bedform's characteristic velocity, sediment, and fluid velocity are differentiated. We developed a new first-order, well-balanced, positivity-preserving, path-preserving, and central wind (WBPP-PCCU) scheme to solve the proposed hyperbolic sediment transport model (HSTM). We used the Averaging Essentially Non-Oscillatory (AENO) reconstruction coupled with the third-order Runge-Kutta Semi-Implicit (SI-RK3) method to achieve second-order accuracy. The balance and positivity of the water depth properties were proven. In this work, a resonance condition is proposed. The model facilitates the application of several other schemes such as Roe, HLLC, HLLEM, PVM (polynomial viscosity matrix), RVM (rational viscosity matrix), which require the diagonalization of the Jacobian matrix. The accuracy, robustness, positivity preservation, and equilibrium properties of the resulting model are evaluated using a series of carefully selected test cases. The proposed model provides an excellent ability to simulate sediment transport in a wide range of coastal environments.

A positivity preserving strategy for entropy stable discontinuous Galerkin discretizations of the compressible Euler and Navier-Stokes equations

High-order entropy-stable discontinuous Galerkin methods for the compressible Euler and Navier-Stokes equations require the positivity of thermodynamic quantities in order to guarantee their well-posedness. In this work, we introduce a positivity limiting strategy for entropy-stable discontinuous Galerkin discretizations constructed by blending high order solutions with a low order positivity-preserving discretization. The proposed low order discretization is semi-discretely entropy stable, and the proposed limiting strategy is positivity preserving for the compressible Euler and Navier-Stokes equations. Numerical experiments confirm the high order accuracy and robustness of the proposed strategy.

Structure-preserving finite volume arbitrary Lagrangian-Eulerian WENO schemes for the shallow water equations

This paper develops the structure-preserving finite volume weighted essentially non-oscillatory (WENO) hybrid schemes for the shallow water equations under the arbitrary Lagrangian-Eulerian (ALE) framework, dubbed as ALE-WENO schemes. The WENO hybrid reconstruction is adopted on moving meshes, which distinguishes the smooth, non-smooth, and transition stencils by a simple smoothness detector. To maintain the positivity preserving and the well-balanced properties of the ALE-WENO schemes, we adapt the positivity preserving limiter and the well-balanced approaches on static meshes to moving meshes. The rigorous theoretical analysis and numerical examples demonstrate the high order accuracy and positivity-preserving property of the schemes under the ALE framework. For the well-balanced schemes, it is successful in the unique exact equilibrium preservation and capturing small perturbations of the hydrostatic state well without numerical oscillations near the discontinuity. Moreover, our ALE-WENO hybrid schemes have an advantage over the simulations on static meshes due to the higher resolution interface tracking of the fluid motion.

A positivity-preserving numerical algorithm for stochastic age-dependent population system with Lévy noise in a polluted environment

This paper develops a model for describing a stochastic age-dependent population system (SADPS) with Lévy noise in a polluted environment. As the model includes a nonlinear drift term and Lévy noise, it is difficult, if not impossible, to derive an analytical solution of the model. This study thus aims at developing a numerical algorithm that is based on a modified truncated Euler-Maruyama (EM) method and positivity preserving. We first study global positivity for the solution of the SADPS model in a polluted environment. Subsequently, we present a semi-discrete Galerkin finite element method in an age variable, and construct a full-discrete numerical approximation using modified truncated EM scheme on time to preserve positivity. Under certain suitable regularity conditions, estimates of convergence error of the full-discrete and a semi-discrete scheme are presented. Finally, a numerical example is given to illustrate the theoretical results.

English

English

Implementation of the HLL-GRP solver for multidimensional ideal MHD simulations based on finite volume method

The HLL solver for the generalized Riemann problem (HLL-GRP solver) is attractive because it is easy to implement for 1D hyperbolic problems and performs well. Based on the rotational invariance of multidimensional magnetohydrodynamics (MHD) equations, this paper provides a means of implementing the 1D HLL-GRP solver into the multidimensional ideal MHD simulations within the framework of the finite volume method. Furthermore, the Godunov-Powell source term and the locally divergence-free WLS-ENO reconstruction are adopted to constrain the divergence of the magnetic field. To maintain positive density and pressure, a positivity-preserving (PP) limiter is performed on the reconstructed polynomials of density and pressure. Some comparisons with other methods are listed in 1D, 2D and 3D MHD benchmark cases. The numerical results demonstrate that the scheme is stable and can well control the divergence of the magnetic field.

Higher-order aerodynamic numerical simulations in compressible RANS framework with inverse-ω scale variable

The use of higher-order methods to robustly compute turbulent flows governed by the multi-equation turbulence model in Reynolds-averaged Navier-Stokes (RANS) framework is a challenging topic in the computational fluid dynamics (CFD) community. In this paper, we propose two aspects of measures involving the physical model and the numerical scheme. Firstly, a λ-scale (=1/ωn) is derived instead of the well-known ω in SSG/LRR Reynolds-stress model (RSM) and SST eddy viscosity model (EVM) due to its natural boundary conditions for viscous walls and help to the positivity preservation of dissipation rate ε in turbulence equations when n=1/(2m),m=1,2,3.... Secondly, a geometric-conservation-preservation negativity-correction (GCPNC) method is proposed on the implicit time integration procedure with right-hand-side (RHS) term discretized by high-order weighted compact nonlinear schemes (WCNSs). It can guarantee the positivity of normal Reynolds stresses or turbulent kinetic energy. Numerical tests with real-life compressible flows are reported to demonstrate the robustness and accuracy of the present methods. Relevant cases include the supersonic square duct flow, the flows over delta wing, wing-body configuration and high-lift configuration.

A simulation expressivity of the quenching phenomenon in a two-sided space-fractional diffusion equation

The aims of this paper are to investigate and propose a numerical approximation for a quenching type diffusion problem associated with a two-sided Riemann-Liouville space-fractional derivative. The approach adopts weighted Grünwald formulas for suitable spatial discretization. An implicit Crank-Nicolson scheme combined with adaptive technology is then implemented for a temporal integration. Monotonicity, positivity preservation and linearized stability are proved under suitable constraints on spatial and temporal discretization parameters. Two specially designed simulation experiments are presented for illustrating and outreaching properties of the numerical method constructed. Connections between the two-sided fractional differential operator and critical values including quenching time, critical length and quenching location are investigated.