Time Step Restriction Research Articles

A fully implicit WENO scheme on stratigraphic and unstructured polyhedral grids

Many authors have used higher-order spatial discretizations to reduce numerical diffusion, which can be particularly pronounced when simulating EOR processes involving active chemical substances that are transported by linear or weakly nonlinear waves. Most high-resolution methods reported in the literature are based on explicit temporal discretizations. This imposes severe time-step restrictions when applied to the type of grids seen in industry-standard simulation models of real assets, which usually have orders-of-magnitude variations in porosities and Darcy velocities that necessitate the use of implicit discretization. Herein, we propose a second-order WENO discretization suitable for complex grids with polyhedral cell geometries, unstructured topologies, large aspect ratios, and large variations in interface areas. The WENO scheme is developed as part of a standard, fully implicit formulation that solves for pressure and transported quantities simultaneously. We investigate the accuracy and utility of the WENO scheme for a series of test cases that involve corner-point and 2D/3D Voronoi grids and black-oil and compositional flow models.

LDG Methods for Reaction-diffusion Systems with Application of Krylov Implicit Integration Factor Methods

In this paper, we present an efficient fully-discrete local discontinuous Galerkin (LDG) method for nonlinear reaction-diffusion systems, which are often used as mathematical models for many physical and biological applications. We can derive numerical approximations not only for solutions but also for their gradients at the same time, while most of methods derive numerical solutions only. And due to the strict time-step restriction ($\Delta t = O(h^2_{\min})$) of explicit schemes for stability, we introduce the implicit integration factor (IIF) method based on Krylov subspace approximation, in which the time step can be taken as $\Delta t = O(h_{\min})$. Moreover, the method allows us to compute element by element and avoid solving a global system of nonlinear algebraic equations as the standard implicit schemes do, which can reduce the computational cost greatly. Numerical experiments about the reaction-diffusion equations with exact solutions and the well-studied Schnakenberg system are conducted to illustrate the accuracy, capability and advantages of the method.

Critical time-step for central difference integration schemes in discrete methods: Translational and rotational degrees of freedom

The explicit central-difference time integration scheme is widely used in discrete methods. However, restrictions on the size of the time step apply. We revisit the problem of a network of particles with translational and rotational degrees of freedom. We apply a discrete Fourier transform to the equations of motion. By studying the eigenvalues of the amplification matrix, we derived a closed form, sharp stability limit that applies to any network. The time-step limit is compared with previous work (Otsubo et al. 2017 and O’Sullivan et al. 2004) for common network configurations. Numerical simulation is used for a certain class of networks in two-dimensional and three-dimensional spaces, and good agreement is observed between the analytical critical time-step and the numerical solutions.

A truly forward semi-Lagrangian WENO scheme for the Vlasov-Poisson system

A class of modular high order forward semi-Lagrangian (FSL) schemes with a number of advantages (outlined below) for linear advection equations is extended to have high resolution, in the sense that sharp gradients can be captured without loss of accuracy, by means of weighted essentially non-oscillatory derivative (DWENO) calculations based on finite differences. Principally, our approach (the convected scheme, CS) differs from conventional high order FSL schemes in the means by which it reaches higher order accuracy. The high order CS progressed in this work reserves use of a low order histogram representation of the solution and a mass-conservative volume-weighted projection operator which ordinarily produces a second order accurate solution. A higher order version is obtained by compensating for the amount by which the CS solution and the exact solution disagree at every order in a Taylor series sense. Instead of including these error terms algebraically, which changes the projection operator itself, the novel idea of the CS is to incorporate this information geometrically as a commensurate adjustment to the volume-weighted proportion assigned to each cell. We leverage DWENO approximations to compute these Taylor series error terms more sensitively, and show the corresponding WENO convected schemes (WENO-CS) achieve the desired dual high order and high resolution behavior. Fifth, seventh, and ninth order WENO-CS base solvers are presented for the 1D constant speed advection equation prototype, and are employed in a Strang splitting approach for the numerical solution of linear and nonlinear two-dimensional hyperbolic equations. Results are showcased for 2D linear transport and rigid body rotation, as well as for classic benchmark problems for the (1+1)-dimensional Vlasov–Poisson system including Landau damping, and the two-stream instability. The schemes enjoy the usual benefits of semi-Lagrangian schemes (e.g. no CFL time step restriction) while further claiming distinct advantages afforded to forward schemes including automatic mass conservation, compact support, and therefore straightforward parallelization.

Staggered-grid entropy-stable multidimensional summation-by-parts discretizations on curvilinear coordinates

The entropy conservative/stable staggered grid tensor-product algorithm of Parsani et al. [1] is extended to multidimensional SBP discretizations. The required SBP preserving interpolation operators are proven to exist under mild restrictions and the resulting algorithm is proven to be entropy conservative/stable as well as elementwise conservative. For 2-dimensional simplex elements, the staggered grid algorithm is shown to be more accurate and have a larger maximum time step restriction as compared to the collocated algorithm. The staggered algorithm significantly reduces the number of (computationally expensive) two-point flux evaluations, which is potentially important for both explicit and implicit time-marching schemes. Furthermore, the staggered algorithm requires fewer degrees of freedom for comparable accuracy, which has favorable implications for implicit time-marching schemes.

An asymptotic preserving scheme for the two-dimensional shallow water equations with Coriolis forces

We consider the two-dimensional Saint-Venant system of shallow water equations with Coriolis forces. We focus on the case of a low Froude number, in which the system is stiff and conventional explicit numerical methods are extremely inefficient and often impractical. Our goal is to design an asymptotic preserving (AP) scheme, which is uniformly asymptotically consistent and stable for a broad range of (low) Froude numbers. The goal is achieved using the flux splitting proposed in [Haack et al., Commun. Comput. Phys., 12 (2012), pp. 955–980] in the context of isentropic Euler and Navier-Stokes equations. We split the flux into the stiff and nonstiff parts and then use an implicit-explicit approach: apply an explicit hyperbolic solver (we use the second-order central-upwind scheme) to the nonstiff part of the system while treating the stiff part of it implicitly. Moreover, the stiff part of the flux is linear and therefore we reduce the implicit stage of the proposed method to solving a Poisson-type elliptic equation, which is discretized using a standard second-order central difference scheme.We conduct a series of numerical experiments, which demonstrate that the developed AP scheme achieves the theoretical second-order rate of convergence and the time-step stability restriction is independent of the Froude number. This makes the proposed AP scheme an efficient and robust alternative to fully explicit numerical methods.

Third order maximum-principle-satisfying DG schemes for convection-diffusion problems with anisotropic diffusivity

For a class of convection-diffusion equations with variable diffusivity, we construct third order accurate discontinuous Galerkin (DG) schemes on both one and two dimensional rectangular meshes. The DG method with an explicit time stepping can well be applied to nonlinear convection–diffusion equations. It is shown that under suitable time step restrictions, the scaling limiter proposed in Liu and Yu (2014) [23] when coupled with the present DG schemes preserves the solution bounds indicated by the initial data, i.e., the maximum principle, while maintaining uniform third order accuracy. These schemes can be extended to rectangular meshes in three dimension. The crucial for all model scenarios is that an effective test set can be identified to verify the desired bounds of numerical solutions. This is achieved mainly by taking advantage of the flexible form of the diffusive flux and the adaptable decomposition of weighted cell averages. Numerical results are presented to validate the numerical methods and demonstrate their effectiveness.

Sharp numerical simulation of incompressible two-phase flows

We present a numerical method for simulating incompressible immiscible fluids, in two and three spatial dimensions. It is constructed as a modified pressure correction projection method on adaptive non-graded Oc/Quadtree Cartesian grids, using the level-set framework to capture the moving interface between the two fluids. The sharp treatment of the interface position, of the fluid parameter discontinuities, and of the interfacial jump conditions ensures convergence in the L∞-norm. Using a novel construction for the pressure guess, we are able to alleviate the standard time step restriction incurred by capillary forces. The solver is validated numerically and employed to simulate the dynamics of physically relevant problems such as rising bubbles and viscous droplets in electric fields.

Unconditional optimal error estimate of the projection/Lagrange-Galerkin finite element method for the Boussinesq equations

This paper provides an unconditional optimal convergence of a fractional-step method for solving the Boussinesq equations. In this method, the convection is treated by the Lagrange-Galerkin technique, whereas the diffusion and the incompressibility are treated by the projection method. There are lots of authors who worked on this method, and some authors gave the error estimate of this method. But, to our best knowledge, the error estimate for this method is under certain time-step restrictions. In this paper, we prove that the methods are stable almost unconditionally, i.e., when τ and h are smaller than a given constant. The basic idea of our analysis is splitting the error function into three terms, one term between the finite element solution and the projection, the other term between the projection and the time-discrete solution, the third term between the time-discrete solution and the exact solution, and giving the error estimates for each term respectively. Then, we obtain the optimal error estimates in L2 and H1-norm for the velocity and L2-norm for the pressure. In order to show the efficiency of our method, some numerical results are presented.

Discontinuous Galerkin discretizations of the Boltzmann–BGK equations for nearly incompressible flows: Semi-analytic time stepping and absorbing boundary layers

We present an efficient nodal discontinuous Galerkin method for approximating nearly incompressible flows using the Boltzmann equations. The equations are discretized with Hermite polynomials in velocity space yielding a first order conservation law. A stabilized unsplit perfectly matching layer (PML) formulation is introduced for the resulting nonlinear flow equations. The proposed PML equations exponentially absorb the difference between the nonlinear fluctuation and the prescribed mean flow. We introduce semi-analytic time discretization methods to improve the time step restrictions in small relaxation times. We also introduce a multirate semi-analytic Adams–Bashforth method which preserves efficiency in stiff regimes. Accuracy and performance of the method are tested using distinct cases including isothermal vortex, flow around square cylinder, and wall mounted square cylinder test cases.

A massively parallel semi-Lagrangian solver for the six-dimensional Vlasov–Poisson equation

This article presents an optimized and scalable semi-Lagrangian solver for the Vlasov–Poisson system in six-dimensional phase space. Grid-based solvers of the Vlasov equation are known to give accurate results. At the same time, these solvers are challenged by the curse of dimensionality resulting in very high memory requirements, and moreover, requiring highly efficient parallelization schemes. In this article, we consider the 6-D Vlasov–Poisson problem discretized by a split-step semi-Lagrangian scheme, using successive 1-D interpolations on 1-D stripes of the 6-D domain. Two parallelization paradigms are compared, a remapping scheme and a domain decomposition approach applied to the full 6-D problem. From numerical experiments, the latter approach is found to be superior in the massively parallel case in various respects. We address the challenge of artificial time step restrictions due to the decomposition of the domain by introducing a blocked one-sided communication scheme for the purely electrostatic case and a rotating mesh for the case with a constant magnetic field. In addition, we propose a pipelining scheme that enables to hide the costs for the halo communication between neighbor processes efficiently behind useful computation. Parallel scalability on up to 65,536 processes is demonstrated for benchmark problems on a supercomputer.

Linearized Galerkin FEMs for Nonlinear Time Fractional Parabolic Problems with Non-smooth Solutions in Time Direction

A Newton linearized Galerkin finite element method is proposed to solve nonlinear time fractional parabolic problems with non-smooth solutions in time direction. Iterative processes or corrected schemes become dispensable by the use of the Newton linearized method and graded meshes in the temporal direction. The optimal error estimate in the $$L^2$$ -norm is obtained without any time step restrictions dependent on the spatial mesh size. Such unconditional convergence results are proved by including the initial time singularity into concern, while previous unconditional convergent results always require continuity and boundedness of the temporal derivative of the exact solution. Numerical experiments are conducted to confirm the theoretical results.

A robust numerical model for shallow water governing solute transport with wet/dry interfaces

In this study, a robust numerical model is developed to guarantee that the C-property for both hydrodynamic and solute transport models and positivity of both water depth and solute concentration. Moreover, special attention is paid to satisfying such properties in the unstructured triangular cells near wet–dry interfaces. In the developed numerical model, the finite-volume central-upwind scheme is adopted to approximate the advective fluxes, and Green’s theorem is applied to estimate the diffusive fluxes. The C-property for hydrodynamic model near wet–dry area is achieved by introducing an auxiliary free surface variable used in the piecewise linear reconstruction and discretizing the bed slope source term in a special form. And the C-property for solute transport is guaranteed by reconstructing secondary variables instead of the primary variables within the frame of central-upwind method, and using the coupling approach for the governing system. The positivity of the solute concentration is satisfied using a proposed diluting time step restriction and the positivity of the water depth is guaranteed by a draining time technique. The desired properties of the proposed numerical model are verified by several numerical experiments.

An explicit–implicit finite element model for the numerical solution of incompressible Navier–Stokes equations on moving grids

In this paper an efficient mesh-moving Finite Element model for the simulation of the incompressible flow problems is proposed. The model is based on a combination of the explicit multi-step scheme (Runge–Kutta) with an implicit treatment of the pressure. The pressure is decoupled from the velocity and is solved for only once per time step minimizing the computational cost of the implicit step. Novel solution algorithm alleviating time step restrictions faced by the majority of the former Lagrangian approaches is presented. The method is examined with respect to its space and time accuracy as well as the computational cost. Two numerical examples are solved: one involving a problem on a domain with fixed boundaries and the other one dealing with a free surface flow. It is shown that the method can be easily parallelized.

Fully implicit interface tracking for a viscous drop under simple shear

In this article we present a novel 3D implicit interface tracking method for sharp interfaces with interfacial tension in the creeping flow regime, employing the finite element method. The interface nodes are allowed to move only in the normal direction and thus remeshing can be avoided, most of the times. The implicit method allows us to overcome certain time step limitations imposed by the mesh capillary time. To validate our method, we use a fairly simple and very well understood problem of a single viscous drop suspended in a viscous matrix that deforms under an applied shear rate. This problem was first studied by Taylor [1] and has been extensively reviewed by Rallison [2] and Stone [3]. The second moment of inertia tensor was used to compute the deformation parameter D and the inclination angle θ and the results are compared to the theory for small deformations of Taylor [1].

Projective and telescopic projective integration for the nonlinear BGK and Boltzmann equations

We present high-order, fully explicit projective integration schemes for nonlinear collisional kinetic equations such as the BGK and Boltzmann equation. The methods first take a few small (inner) steps with a simple, explicit method (such as direct forward Euler) to damp out the stiff components of the solution. Then, the time derivative is estimated and used in an (outer) Runge-Kutta method of arbitrary order. The procedure can be recursively repeated on a hierarchy of projective levels to construct telescopic projective integration methods. Based on the spectrum of the linearized collision operator, we deduce that the computational cost of the method is essentially independent of the stiffness of the problem: with an appropriate choice of inner step size, the time step restriction on the outer time step, as well as the number of inner time steps, is independent of the stiffness of the (collisional) source term. In some cases, the number of levels in the telescopic hierarchy depends logarithmically on the stiffness. We illustrate the method with numerical results in one and two spatial dimensions.

A parabolic level set reinitialisation method using a discontinuous Galerkin discretisation

Level set reinitialisation is a part of the level set methodology which allows one to generate, at any point during level set evolution, a level set function which is a signed distance function to its own zero isocontour. Whilst not in general a required condition, maintaining the level set function as a signed distance function is often desirable as it removes a known source of numerical instability. This paper presents a novel level set reinitialisation method based on the solution of a nonlinear parabolic PDE. The PDE is discretised using a symmetric interior penalty discontinuous Galerkin method in space, and an implicit Euler method in time. Also explored are explicit and semi-implicit time discretisations, however, numerical experiments demonstrate that such methods suffer from severe time step restrictions, leading to prohibitively large numbers of iterations required to achieve convergence. The proposed method is shown to be high-order accurate through a number of numerical examples. More specifically, the presented experimental orders of convergence align with the well established optimal convergence rates for the symmetric interior penalty method; that is the error in the L2 norm decreases proportionally to hp+1 and the error in the DG norm decreases proportionally to hp.

Improvements in highly viscous fluid simulation using a fully implicit SPH method

This study describes the application of two main improvements in highly viscous fluid simulations using the smoothed particle hydrodynamics (SPH) method: an implicit time integration scheme to overcome the problem of impractically small time step restriction and the introduction of air ghost particles to fix problems regarding the free surface treatment. This study adopts the incompressible SPH as a basis for the implementation of these improvements, which guarantees a stable and accurate pressure distribution. We verified the proposed implicit time integration scheme with simulations of pipe flow and the free surface treatment with a simple hydrostatic problem. As a result, the free surface of the hydrostatic problem became very smooth and stable. In addition, we conducted a variety of dam-break simulations to validate this proposed SPH method, as well as to analyze the density and divergence error. Finally, we demonstrate the potential of this method with the highly viscous vertical jet flow over a horizontal plate test, which features a complex viscous coiling behavior.

Fully implicit time integration in truly incompressible SPH

In the last years, the Smoothed Particle Hydrodynamics (SPH) method was introduced in new fields of engineering applications where highly viscous flow like polymer flow or very small geometries, like in micro-flow, are involved. In most of these applications it was not possible to use realistic fluid parameters, larger simulation domains or higher resolution because of restriction of time step due to the viscous time step criterion of the explicit time integration scheme. The computational effort was too high even for highly scalable SPH codes. In this article, we present a first-order implicit time stepping scheme for truly incompressible SPH (ISPH) to eliminate the viscous time step criterion. We propose a consistent time stepping approach where both velocity and particle position are solved implicitly and compare the results to traditional semi-implicit ISPH and implicit schemes where only velocity is integrated implicitly. We study Poiseuille flow, Taylor–Green flow and Rayleigh–Taylor instability and compare results of the implicit schemes to ISPH. Energy conservation is investigated using Taylor–Green vortex to highlight differences between the approaches. We find that the fully implicit time integration scheme is numerically stable. Since no error estimation of accuracy is used, the larger time step size leads to deviations of the trajectories. In future work, higher-order time integration schemes as well as error estimators should be investigated.

An implicit unified gas-kinetic scheme for unsteady flow in all Knudsen regimes

The unified gas-kinetic scheme (UGKS) is a direct modeling method for multiple scale transport. Based on the ratio of time step to the particle collision time, the local evolution solution on the mesh size and time step scales is used in the construction of the multiscale method. For a flow problem covering multiple flow regimes, such as the hypersonic flow around a flying vehicle in near space, the UGKS is able to capture the highly compressed Navier-Stokes solution in one region and fully expanded free molecular flow in another region, with significant variations of the ratio between the time step and the local particle collision time around the vehicle. For an explicit UGKS, the time step in the whole computational domain is determined by the CFL condition. With implicit and multigrid techniques, the efficiency of the UGKS [1,2] has been improved by two orders of magnitude for steady state computation. However, for unsteady flow computation, due to the CFL condition the global time step used in the explicit UGKS may be limited by the smallest cell size in the computational domain. As a result, for a largely stretched non-uniform mesh the global time step becomes very small and the ratio of the time step to the local particle collision time may get a very small value. Under such a circumstance, even though the UGKS is a multiscale method, the real physics represented in the explicit UGKS may be constrained to the kinetic scale transport only, and the advantage of the multiscale nature in UGKS has not been fully utilized. In order to solve the multiscale unsteady flow problem efficiently, the time step restriction from a global CFL condition has to be released. In this paper, we will develop an implicit UGKS (IUGKS) for unsteady flows by alternatively solving the macroscopic and microscopic governing equations within a time step iteratively. With a pre-defined uniform large evolution time step, the local CFL number varies greatly in different region, such as on the order 1 in the large numerical cell size region, and 100 in the small cell size region. In order to preserve coherent flow evolution and keep the multiscale nature, the time averaged numerical flux across a cell interface is still evaluated by the explicit UGKS under the local CFL condition. Therefore, the multiscale property of the UGKS modeling has been kept over non-uniform meshes. With improved temporal discretization, the current IUGKS can automatically go back to the explicit UGKS and obtain identical solutions when the time step of the implicit scheme gets to that of an explicit one. Many numerical examples are included to validate the scheme for both continuum and rarefied flows with a large variation of artificially generated mesh size. The IUGKS has a second order accuracy and presents reasonably good results for unsteady flow computation, and its efficiency has been improved by dozens of times in comparison with the explicit UGKS.