Search

In this paper, we propose and analyze a new class of schemes based on a variant of the scalar auxiliary variable (SAV) approaches for gradient flows. Precisely, we construct more robust first and second order unconditionally stable schemes by introducing a new defined auxiliary variable to deal with nonlinear terms in gradient flows. The new approach consists in splitting the gradient flow into decoupled linear systems with constant coefficients, which can be solved using existing fast solvers for the Poisson equation. This approach can be regarded as an extension of the SAV method; see, e.g., Shen et al. (2018) [21] , in the sense that the new approach comes to be the conventional SAV method when α = 0 and removes the boundedness assumption on ∫ Ω F ( ϕ ) d x required by the SAV. The new approach only requires that the total free energy or a part of it is bounded from below, which is more realistic in physically meaningful models. The unconditional stability is established, showing that the efficiency of the new approach is less restricted to particular forms of the nonlinear terms. A series of numerical experiments is carried out to verify the theoretical claims and illustrate the efficiency of our method.

The Keller-Segel equations are widely used for describing chemotaxis in biology. Recently, first-order and second-order approximations for a class of Keller-Segel equations based on the gradient flow structure were proposed in [48] . Mass conservation, positivity and energy stability were proved for the first-order scheme, whereas for the second-order scheme the energy stability was not provided. Besides, an explicit-implicit treatment is performed to a non-convex and non-concave term − χ ρ ϕ , making their decoupled system could only be solved in sequence. In this paper, we propose new BDF schemes of first-order (BDF1) and second-order accuracy (BDF2 and EsBDF2): the coupled term − χ ρ ϕ involved in two equations of ρ and ϕ is fully explicitly treated, thus the discrete schemes could be computed in parallel. For the first-order scheme (BDF1) and the second-order scheme (BDF2), the standard backward differentiation formula is applied to approximate the origin continuous equations with a regularization term added to guarantee the unconditional energy stability. For the BDF1 scheme, the discrete system is proved to be energy stable with a constant restriction for the stabilized parameter. For the EsBDF2 scheme, which is different with the standard BDF2 scheme, the differential operators are carefully handled and some regularization terms are added to provide the energy stability. Several numerical examples are presented to verify the theoretical results.

We propose a Convex Splitting Runge–Kutta (CSRK) scheme which provides a simple unified framework to solve a gradient flow in an unconditionally gradient stable manner. The key feature of the scheme is a combination of a convex splitting method and a specially designed multi-stage two-additive Runge–Kutta method. Our methods are high order accurate in time and assure the gradient (energy) stability for any time step size. We provide detailed proof of the unconditional energy stability and present issues on the practical implementations. We demonstrate the accuracy and stability of the proposed methods using numerical experiments of the Cahn–Hilliard equation.

Multiphase flow in porous media occurs in several disciplines including petroleum reservoir engineering, petroleum systems' analysis, and CO$_2$ sequestration. While simulations often use a fully implicit discretization to increase the time step size, restrictions on the time step often exist due to non-convergence of the nonlinear solver (e.g. Newton's method). Here this problem is addressed for the Buckley-Leverett equations, which model incompressible, immiscible, two-phase flow with no capillary potential. The equations are recast as a gradient flow using the phase-field method, and a convex energy splitting scheme is applied to enable large timesteps, even for high degrees of heterogeneity in permeability and viscosity. By using the phase-field formulation as a homotopy map, the underlying hyperbolic flow equations can be solved with large timesteps. For a heterogeneous test problem, the new homotopy method allows the timestep to be increased by more than six orders of magnitude relative to the unmodified equations while maintaining convergence., Comment: 12 pages, 7 figures, 2 tables

We present a paradigm for developing arbitrarily high order, linear, unconditionally energy stable numerical algorithms for general gradient flow models. We apply the energy quadratization (EQ) technique to reformulate the gradient flow model into an equivalent one with a quadratic free energy and a modified mobility. For any positive integer k > 0 and t n = n Δ t , where Δt is the time step size, we linearize the EQ-reformulated model in ( t n , t n + 1 ] through an extrapolation from numerical solutions already obtained in [ 0 , t n ] so that the approximation error is in the order of O ( Δ t k ) . Then we employ an s-stage algebraically stable Runge-Kutta method to discretize the linearized model in ( t n , t n + 1 ] . For the spatial discretization, we use the Fourier pseudo-spectral method to match the order of accuracy in time. The resulting fully discrete scheme is linear, unconditionally energy stable, uniquely solvable, and can reach arbitrarily high order. Furthermore, we present a family of linear schemes based on prediction-correction methods to complement the new linear schemes. Some benchmark numerical examples are given to demonstrate the accuracy and efficiency of the schemes.

Phase field models are widely used to describe the two-phase system. The evolution of the phase field variables is usually driven by the gradient flow of a total free energy functional. The generalization of the approach to an N phase ( N ź 3 ) system requires some extra consistency conditions on the free energy functional in order for the model to give physically relevant results. A projection approach is proposed for the derivation of a consistent free energy functional for the three-phase Cahn-Hilliard equations. The system is then coupled with the Navier-Stokes equations to describe the three-phase flow on solid surfaces with moving contact line. An energy stable scheme is developed for the three-phase flow system. The discrete energy law of the numerical scheme is proved which ensures the stability of the scheme. We also show some numerical results for the dynamics of triple junctions and four phase contact lines.

In this paper, we present two efficient and spectrally accurate numerical methods for computing the ground and first excited states in Bose-Einstein condensates (BECs). We begin with a review on the gradient flow with discrete normalization (GFDN) for computing stationary states of a nonconvex minimization problem and show how to choose initial data effectively for the GFDN. For discretizing the gradient flow, we use sine-pseudospectral method for spatial derivatives and either backward Euler scheme (BESP) or backward/forward Euler schemes for linear/nonlinear terms (BFSP) for temporal derivatives. Both BESP and BFSP are spectral order accurate for computing the ground and first excited states in BEC. Of course, they have their own advantages: (i) for linear case, BESP is energy diminishing for any time step size where BFSP is energy diminishing under a constraint on the time step size; (ii) at every time step, the linear system in BFSP can be solved directly via fast sine transform (FST) and thus it is extremely efficient, and in BESP it needs to be solved iteratively via FST by introducing a stabilization term and thus it could be efficient too. Comparisons between BESP and BFSP as well as other existing numerical methods are reported in terms of accuracy and total computational time. Our numerical results show that both BESP and BFSP are much more accurate and efficient than those existing numerical methods in the literature. Finally our new numerical methods are applied to compute the ground and first excited states in BEC in one dimension (1D), 2D and 3D with a combined harmonic and optical lattice potential for demonstrating their efficiency and high resolution.

A noniterative implicit space-marching method based on an adaptive grid approach is developed for solving the 2D steady state Euler equations describing supersonic gas flows without streamwise separation. A grid adapted to the solution is generated by using a variational method optimizing such important properties of a grid as smoothness, orthogonality, and grid cell volume variation simultaneously. The weight function in the integral of adaptation is constructed, so that the grid lines are accumulated near strong gradients and curvatures of the solution curve. A special treatment of the boundary points is proposed to determine the location of nodes near slope discontinuities of the boundary mesh line. The noniterative implicit space-marching method of Y. C. Vigneronet al. is applied to solve the Euler equations written in an arbitrary curvilinear system of coordinates. The streamwise and transverse derivatives in the governing equations are approximated by the second-order upwind Richardson and second-order symmetric TVD schemes, respectively. Implicit boundary condition procedure based on the theory of characteristics for hyperbolic systems of equations is employed. Numerical calculations show that the resolution of strong gradient flow fields can significantly be improved by using the grid adapted to the solution, while the number of grid nodes is the same.

A coupled level set method for the motion of multiple junctions (of, e.g., solid, liquid, and grain boundaries), which follows the gradient flow for an energy functional consisting of surface tension (proportional to length) and bulk energies (proportional to area), is developed. The approach combines the level set method of S. Osher and J. A. Sethian with a theoretical variational formulation of the motion by F. Reitich and H. M. Soner. The resulting method uses as many level set functions as there are regions and the energy functional is evaluated entirely in terms of level set functions. The gradient projection method leads to a coupled system of perturbed (by curvature terms) Hamilton?Jacobi equations. The coupling is enforced using a single Lagrange multiplier associated with a constraint which essentially prevents (a) regions from overlapping and (b) the development of a vacuum. The numerical implementation is relatively simple and the results agree with (and go beyond) the theory as given in 12. Other applications of this methodology, including the decomposition of a domain into subregions with minimal interface length, are discussed. Finally, some new techniques and results in level set methodology are presented.

A suite of seven test cases is proposed for the evaluation of numerical methods intended for the solution of the shallow water equations in spherical geometry. The shallow water equations exhibit the major difficulties associated with the horizontal dynamical aspects of atmospheric modeling on the spherical earth. These cases are designed for use in the evaluation of numerical methods proposed for climate modeling and to identify the potential trade-offs which must always be made in numerical modeling. Before a proposed scheme is applied to a full baroclinic atmospheric model it must perform well on these problems in comparison with other currently accepted numerical methods. The cases are presented in order of complexity. They consist of advection across the poles, steady state geostrophically balanced flow of both global and local scales, forced nonlinear advection of an isolated low, zonal flow impinging on an isolated mountain, Rossby-Haurwitz waves, and observed atmospheric states. One of the cases is also identified as a computer performance/algorithm efficiency benchmark for assessing the performance of algorithms adapted to massively parallel computers.

We consider a class of time-dependent second order partial differential equations governed by a decaying entropy. The solution usually corresponds to a density distribution, hence positivity (non-negativity) is expected. This class of problems covers important cases such as Fokker–Planck type equations and aggregation models, which have been studied intensively in the past decades. In this paper, we design a high order discontinuous Galerkin method for such problems. If the interaction potential is not involved, or the interaction is defined by a smooth kernel, our semi-discrete scheme admits an entropy inequality on the discrete level. Furthermore, by applying the positivity-preserving limiter, our fully discretized scheme produces non-negative solutions for all cases under a time step constraint. Our method also applies to two dimensional problems on Cartesian meshes. Numerical examples are given to confirm the high order accuracy for smooth test cases and to demonstrate the effectiveness for preserving long time asymptotics.