Semi-implicit Time Discretization Research Articles

A new class of efficient high order semi-Lagrangian IMEX discontinuous Galerkin methods on staggered unstructured meshes

In this paper we present a new high order semi-implicit discontinuous-Galerkin (DG) scheme on two-dimensional staggered triangular meshes applied to different nonlinear systems of hyperbolic conservation laws such as advection-diffusion models, incompressible Navier-Stokes equations and natural convection problems. While the temperature and pressure field are defined on a triangular main grid, the velocity field is defined on a quadrilateral edge-based staggered mesh. A semi-implicit time discretization is proposed, which separates slow and fast time scales by treating them explicitly and implicitly, respectively. The nonlinear convection terms are evolved explicitly using a semi-Lagrangian approach, whereas we consider an implicit discretization for the diffusion terms and the pressure contribution. High-order of accuracy in time is achieved using a new flexible and general framework of IMplicit-EXplicit (IMEX) Runge-Kutta schemes specifically designed to operate with semi-Lagrangian methods. To improve the efficiency in the computation of the DG divergence operator and the mass matrix, we propose to approximate the numerical solution with a less regular polynomial space on the edge-based mesh, which is defined on two sub-triangles that split the staggered quadrilateral elements. This allows for a fast computation of the DG operators and matrices without any need to store them for each mesh element. Due to the implicit treatment of the fast scale terms, the resulting numerical scheme is unconditionally stable for the considered governing equations because the semi-Lagrangian approach is the only explicit discretization for convection terms that does not need any stability restriction on the maximum admissible time step. Contrarily to a genuinely space-time discontinuous-Galerkin scheme, the IMEX discretization permits to preserve the symmetry and the positive semi-definiteness of the arising linear system for the pressure that can be solved at the aid of an efficient matrix-free implementation of the conjugate gradient method. We present several convergence results, including nonlinear transport and density currents, up to third order of accuracy in both space and time.

A mathematical analysis to the approximate weak solution of the Smagorinsky Model for different flow regimes

This study delves into the numerical approximation of non-stationary Navier-Stokes equations within turbulent regimes, employing the Smagorinsky Model (SM). By treating the model as inherently discrete, we implement a semi-implicit time discretization using the Euler method. This approach includes comprehensive stability analyses, applicable to a spectrum of flow regimes, and an exploration of the asymptotic energy balance dynamics during fluid movements. The primary contribution of this study is found in its methodical approach to the numerical approximation of non-stationary Navier-Stokes equations within turbulent regimes using the Smagorinsky Model (SM). The adoption of a semi-implicit time discretization with the Euler method, coupled with a meticulous analysis of energy balance, establishes a robust foundation adaptable to diverse flow conditions.

Boundary effects on wave trains in the Exner model of sedimental transport

In this work we compute the numerical solution of the Exner model of sedimentation when a train of waves is imposed at the inflow boundary (E. Macca in Shock-Capturing methods: Well-Balanced Approximate Taylor and Semi-Implicit schemes. PhD thesis, Universit‘a degli Studi di Palermo, Palermo, (2022) and Major Revision, (2023)). The numerical solver is a second order finite-volume scheme, with semi-implicit time discretization based on Implicit-Explicit (IMEX) schemes, which guarantees better stability properties than explicit ones, still at a lower cost than fully implicit schemes. We show the effect of spurious reflected waves generated at the outflow edge of the computational domain, propose two remedies, and show how such spurious effects can be reduced by suitable non-reflecting boundary conditions.

A positivity-preserving implicit-explicit scheme with high order polynomial basis for compressible Navier–Stokes equations

In this paper, we are interested in constructing a scheme solving compressible Navier–Stokes equations, with desired properties including high order spatial accuracy, conservation, and positivity-preserving of density and internal energy under a standard hyperbolic type CFL constraint on the time step size, e.g., Δt=O(Δx). Strang splitting is used to approximate convection and diffusion operators separately. For the convection part, i.e., the compressible Euler equation, the high order accurate postivity-preserving Runge–Kutta discontinuous Galerkin method can be used. For the diffusion part, the equation of internal energy instead of the total energy is considered, and a first order semi-implicit time discretization is used for the ease of achieving positivity. A suitable interior penalty discontinuous Galerkin method for the stress tensor can ensure the conservation of momentum and total energy for any high order polynomial basis. In particular, positivity can be proven with Δt=O(Δx) if the Laplacian operator of internal energy is approximated by the Qk spectral element method with k=1,2,3. So the full scheme with Qk (k=1,2,3) basis is conservative and positivity-preserving with Δt=O(Δx), which is robust for demanding problems such as solutions with low density and low pressure induced by high-speed shock diffraction. Even though the full scheme is only first order accurate in time, numerical tests indicate that higher order polynomial basis produces much better numerical solutions, e.g., better resolution for capturing the roll-ups during shock reflection.

Analysis of dissipation operators that damp spurious modes while maintaining discrete approximate geostrophic equilibriums for the B-grid staggered scheme on triangular meshes

We consider the lowest order - so called B-grid - discretization on triangular cells which is common in the ocean or atmosphere simulation science, with pressure field at the vertices of the triangles and the velocity field at the cell centers. It is known that this discretization presents spurious modes due to a too large number of velocities unknowns and that these modes should be damped; since preservation of approximate geostrophic equilibrium (GE) may be of primary importance in the simulations, we propose and analyze damping operators that, when applied to the linearized constant coefficient configuration, exactly preserve discrete versions of the GEs. We also analyze their effects in more general configurations and their possible coupling with more standard damping operators. These strategies cover triangular variants of the Low Froude (LF) scheme [1], the Apparent Topography (AT) scheme [2,3], as well as a new modification of it. The LF scheme has no numerical diffusion terms in the pressure equation and damps the normal velocity jumps between adjacent cells in the momentum equation. In the AT scheme, numerical diffusion is added in the pressure equation in such a way that it does not impact the discrete GE. However, since the AT scheme alone cannot be proved to be damping through energy estimates and does not preserve the space which is orthogonal to the discrete GE, we also study what we call the Modified Apparent Topography (MAT) that satisfies these important properties. We also suggest semi-implicit time discretizations that preserve these properties. Extensions to the full non-linear equations are also provided. Numerical experiments illustrate the performance of the new schemes.

A conservative Galerkin solver for the quasilinear diffusion model in magnetized plasmas

The quasilinear theory describes the resonant interaction between particles and waves with two coupled equations: one for the evolution of the particle probability density function (pdf), the other for the wave spectral energy density (sed). In this paper, we propose a conservative Galerkin scheme for the quasilinear model in three-dimensional momentum space and three-dimensional spectral space, with cylindrical symmetry.We construct an unconditionally conservative weak form, and propose a discretization that preserves the unconditional conservation property, by “unconditional” we mean that conservation is independent of the singular transition probability. The discrete operators, combined with a consistent quadrature rule, will preserve all the conservation laws rigorously. The technique we propose is quite general: it works for both relativistic and non-relativistic systems, for both magnetized and unmagnetized plasmas, and even for problems with time-dependent dispersion relations.We represent the particle pdf by continuous basis functions, and use discontinuous basis functions for the wave sed, thus enabling the application of a positivity-preserving technique. The marching simplex algorithm, which was initially designed for computer graphics, is adopted for numerical integration on the resonance manifold. We introduce a semi-implicit time discretization, and discuss the stability condition. In addition, we present numerical examples with a “bump on tail” initial configuration, showing that the particle-wave interaction results in a strong anisotropic diffusion effect on the particle pdf.

Abstract nonlinear evolution inclusions of second order with applications in visco-elasto-plasticity

Existence of strong solutions of an abstract Cauchy problem for a class of doubly nonlinear evolution inclusion of second order is established via a semi-implicit time discretization method. The principal parts of the operators acting on u and u′ are multi-valued subdifferential operators and are discretized implicitly. A non-variational and non-monotone perturbation acting nonlinearly on u and u′ is allowed and discretized explicitly in time. The convergence of a variational approximation scheme is established using methods from convex analysis. In addition, it is proven that the solution satisfies an energy-dissipation equality. Applications of the abstract theory to various examples, e.g., a model in visco-elastic-plasticity, are provided.

Convergence of a scheme for an elastic flow with tangential mesh movement

Elastic flow for closed curves can involve significant deformations. Mesh-based approximation schemes require tangentially redistributing vertices for long-time computations. We present and analyze a method that uses the Dirichlet energy for this purpose. The approach effectively also penalizes the length of the curve, and equilibrium shapes are equivalent to stationary points of the elastic energy augmented with the length functional. Our numerical method is based on linear parametric finite elements. Following the lines of Deckelnick and Dziuk [Math. Comp. 78 (2009) 645–671] we prove convergence and establish error estimates, noting that the addition of the Dirichlet energy simplifies the analysis in comparison with the length functional. We also present a simple semi-implicit time discretization and discuss some numerical results that support the theory.

Predicting the effect of a combination drug therapy on the prostate tumor growth via an improvement of a direct radial basis function partition of unity technique for a diffuse-interface model

Chemotherapy is usually applied to treat advanced prostate cancer that cancer cells spread outside the prostate gland. The treatment uses cytotoxic drugs to target cells that grow and divide quickly. On the other hand, the growth of such cancerous tumors depends on angiogenesis. In this paper, we numerically study a diffuse-interface model in a two-dimensional space related to the therapies of prostate cancer. The proposed model describes the tumor growth driven by a generic nutrient and producing the prostate-specific antigen. More precisely, the effect of cytotoxic chemotherapy in the model is evaluated by considering a time-dependent function in the tumor dynamics. Also, another function related to the antiangiogenic therapy is considered to show the reducing intratumoral nutrient supply in the nutrient dynamics. Here, a meshless approximation, i.e., a generalized form of the direct radial basis function partition of unity (D-RBF-PU) method is presented for finding the numerical simulations of this model utilizing in medical oncology. The method uses the lower number of trial points in each patch than the original D-RBF-PU scheme for approximating the trial function per test point. Hence, the time complexity of the method is less than the D-RBF-PU technique. Besides, a semi-implicit time discretization of order 1 has been used to deal with the time variable. Consequently, a linear system of algebraic equations could be solved iteratively per time step by the use of the biconjugate gradient stabilized method with zero-fill incomplete lower–upper preconditioner. Finally, the obtained results without using any adaptive algorithm demonstrate the response of the prostate tumor growth to the chemotherapy, antiangiogenic therapy and a combined therapy.

A projection-based, semi-implicit time-stepping approach for the Cahn-Hilliard Navier-Stokes equations on adaptive octree meshes

The Cahn-Hilliard Navier-Stokes (CHNS) system provides a computationally tractable model that can be used to effectively capture interfacial dynamics in two-phase fluid flows. In this work we present a semi-implicit, projection-based finite element framework for solving the CHNS system. We use a projection-based semi-implicit time discretization for the Navier-Stokes equation and a fully-implicit time discretization for the Cahn-Hilliard equation. We use a conforming continuous Galerkin (cG) finite element method in space equipped with a residual-based variational multiscale (RBVMS) formulation. Pressure is decoupled using a projection step, which results in two linear positive semi-definite systems for velocity and pressure, instead of the saddle point system of a pressure-stabilized method. All the linear systems are solved using an efficient and scalable algebraic multigrid (AMG) method. We deploy this approach on a massively parallel numerical implementation using parallel octree-based adaptive meshes. The overall approach allows the use of relatively large time steps with much faster time-to-solve than similar fully-implicit methods. We present comprehensive numerical experiments showing detailed comparisons with results from the literature for canonical cases, including the single bubble rise and Rayleigh-Taylor instability.

Mathematical Modeling of Wax Deposition in Field-Scale Crude Oil Pipeline Systems

The formation of solid wax crystals, which interlock and form a gel-like layer on the inner wall of the crude oil pipeline, influences the transportation of waxy crude oil. The deposited layer grows continuously and hardens during the oil transportation, reducing the effective inside diameter of the crude oil pipeline and the flow rate. In extreme cases, the deposited layer may block the crude oil pipeline leading to a loss of production and capital investment. In this paper, wax deposition from multiphase flow in field-scale oil pipeline transport systems has been studied. The novelty of this work is to develop a mathematical model that incorporates water-in-oil emulsions, wax precipitation kinetics, molecular diffusion, and shear dispersion to enable accurate predictions of both the wax deposit growth rate and aging of the deposit. The coupled nonlinear partial differential equations governing the flow are discretized in time by a second-order semi-implicit time discretization scheme based on the Adams-Bashforth and Crank-Nicolson methods, which completely decouples the computation of the governing equations. The resulting temporal schemes are discretized in space by the bivariate spectral collocation method based on Chebyshev-Gauss-Lobatto grid points and simulated in MATLAB software to obtain the profiles of the flow variables. The simulation results are presented in graphical and in tabular forms and discussed. This study found that the deposit thickness is directly proportional to the Reynolds number and inversely proportional to the mass Grashof number, Schmidt number, and Weber number. Deposit aging is rampant during the early stages of wax deposition, after which it stabilizes at a specific value as time elapses. A deposition model to predict the wax deposit thickness and aging is proposed in this study. The findings of this study will help in making informed decisions on the planning of pigging operations, thermal insulation, and other remediation techniques to be applied in controlling wax deposition in field-scale crude oil pipeline systems.

A proximal-gradient algorithm for crystal surface evolution

As a counterpoint to recent numerical methods for crystal surface evolution, which agree well with microscopic dynamics but suffer from significant stiffness that prevents simulation on fine spatial grids, we develop a new numerical method based on the macroscopic partial differential equation, leveraging its formal structure as the gradient flow of the total variation energy, with respect to a weighted \(H^{-1}\) norm. This gradient flow structure relates to several metric space gradient flows of recent interest, including 2-Wasserstein flows and their generalizations to nonlinear mobilities. We develop a novel semi-implicit time discretization of the gradient flow, inspired by the classical minimizing movements scheme (known as the JKO scheme in the 2-Wasserstein case). We then use a primal dual hybrid gradient (PDHG) method to compute each element of the semi-implicit scheme. In one dimension, we prove convergence of the PDHG method to the semi-implicit scheme, under general integrability assumptions on the mobility and its reciprocal. Finally, by taking finite difference approximations of our PDHG method, we arrive at a fully discrete numerical algorithm, with iterations that converge at a rate independent of the spatial discretization: in particular, the convergence properties do not deteriorate as we refine our spatial grid. We close with several numerical examples illustrating the properties of our method, including facet formation at local maxima, pinning at local minima, and convergence as the spatial and temporal discretizations are refined.

High order semi-implicit schemes for viscous compressible flows in 3D

In this article we present a high order cell-centered numerical scheme in space and time for the solution of the compressible Navier-Stokes equations. To deal with multiscale phenomena induced by the different speeds of acoustic and material waves, we propose a semi-implicit time discretization which allows the CFL-stability condition to be independent of the fast sound speed, hence improving the efficiency of the solver. This is particularly well suited for applications in the low Mach regime with a rather small fluid velocity, where the governing equations tend to the incompressible model. The momentum equation is inserted into the energy equation yielding an elliptic equation on the pressure. The class of implicit-explicit (IMEX) time integrators is then used to ensure asymptotic preserving properties of the numerical method and to improve time accuracy. High order in space is achieved relying on implicit finite difference and explicit CWENO reconstruction operators, that ultimately lead to a fully quadrature-free scheme. To relax the severe parabolic restriction on the maximum admissible time step related to viscous contributions, a novel implicit discretization of the diffusive terms is designed. A variational approach based on the discontinuous Galerkin (DG) spatial discretization is devised in order to obtain a discrete cell-centered Laplace operator. High order corner gradients of the velocity field are employed in 3D to derive the Laplace discretization, and the resulting viscous system is proven to be symmetric and positive definite. As such, it can be conveniently solved at the aid of the conjugate gradient method. Numerical results confirm the accuracy and the robustness of the novel schemes in the challenging stiff limit of the governing equations characterized by low Mach numbers.

Mean curvature motion of point cloud varifolds

This paper investigates a discretization scheme for mean curvature motion on point cloud varifolds with particular emphasis on singular evolutions. To define the varifold a local covariance analysis is applied to compute an approximate tangent plane for the points in the cloud. The core ingredient of the mean curvature motion model is the regularization of the first variation of the varifoldviaconvolution with kernels with small stencil. Consistency with the evolution velocity for a smooth surface is proven provided that a sufficiently small stencil and a regular sampling are considered. Furthermore, an implicit and a semi-implicit time discretization are derived. The implicit scheme comes with discrete barrier properties known for the smooth, continuous evolution, whereas the semi-implicit still ensures in all our numerical experiments very good approximation properties while being easy to implement. It is shown that the proposed method is robust with respect to noise and recovers the evolution of smooth curves as well as the formation of singularities such as triple points in 2D or minimal cones in 3D.

Adaptive local discontinuous Galerkin methods with semi-implicit time discretizations for the Navier-Stokes equations

In this paper, we present a mesh adaptation algorithm for the unsteady compressible Navier-Stokes equations under the framework of local discontinuous Galerkin methods coupled with implicit-explicit Runge-Kutta or spectral deferred correction time discretization methods. In both of the two high order semi-implicit time integration methods, the convective flux is treated explicitly and the viscous and heat fluxes are treated implicitly. The remarkable benefits of such semi-implicit temporal discretizations are that they can not only overcome the stringent time step restriction compared with time explicit methods, but also avoid the construction of the large Jacobian matrix as is done for fully implicit methods, thus are relatively easy to implement. To save computing time as well as capture the flow structures of interest accurately, a local mesh refinement (h-adaptive) technique, in which we present detailed criteria for selecting candidate elements and complete strategies to refine and coarsen them, is also applied for the Navier-Stokes equations. Numerical experiments are provided to illustrate the high order accuracy, efficiency and capabilities of the semi-implicit schemes in combination with adaptive local discontinuous Galerkin methods for the Navier-Stokes equations.

Solving phase-field models in the tensor train format to generate microstructures of bicontinuous composites

Phase-field models permit accounting for the underlying physics of the microstructure evolution process when simulating the emergent microstructure of a variety of multiscale materials. As an inherent characteristic, the individual phases (or states) of the material are not known in advance and only emerge upon evolution. To account for this matter, structured grids are typically used when solving phase-field models. In particular for large three-dimensional grids and long time scales, the used memory and computational time can be considerable, instigating research on alternative approaches.The work at hand studies the Tensor Train (TT) format as a possible remedy for phase-field simulations. More precisely, the TT format is a specific tensor format which permits – small so-called TT rank provided – to store huge arrays of data and comes with specific algorithms replacing the linear algebra operations familiar for fully stored arrays.We investigate the Cahn-Hilliard equation and use a semi-implicit discretization in time which is based on solving a single linear system per time step. We compare this strategy with a classic implicit Euler discretization using a Newton scheme as a nonlinear solver. Upon a finite-difference discretization in space, the resulting equations are recast in the Quantics Tensor Train (QTT) format and solved by dedicated linear solvers. We study the performance of the algorithms and investigate the effective properties of the resulting bicontinuous composites via computational homogenization. We show that the Cahn-Hilliard equation can be solved using the TT format and that runtimes scale well with the grid size.

High Order Semi-implicit WENO Schemes for All-Mach Full Euler System of Gas Dynamics

In this paper, we propose high order semi-implicit schemes for the all Mach full Euler equations of gas dynamics. Material waves are treated explicitly, while acoustic waves are treated implicitly, thus avoiding severe CFL restrictions for low Mach flows. High order accuracy in time is obtained by semi-implicit temporal integrator based on the IMEX Runge-Kutta (IMEX-RK) framework. High order in space is achieved by finite difference WENO schemes with characteristic-wise reconstructions adapted to the semi-implicit IMEX-RK time discretization. The schemes are proven to be asymptotic preserving and asymptotically accurate as the Mach number vanishes. Besides, they can well capture discontinuous solutions in the compressible regime, especially for two dimensional Riemann problems. Numerical tests in one and two space dimensions will illustrate the effectiveness of the proposed schemes.

A Stable and High Accurate Numerical Simulator for Convection-Diffusion Equation

This study presents a simulator to obtain numerical solution of convection-diffusion equation. It includes explicit, fully implicit and semi-implicit time discretization techniques. Although the explicit method is efficient for some simple conduction problems, it needs stability criteria. Otherwise, it is unstable especially for high Courant number and large time step size. When using the explicit method, the time step size is chosen with care to get well-posed numerical solution. This requirement sets a serious constraint for the explicit method because choosing small time step size causes the simulation time to become quite long. Even though the implicit method is stable for large time step size and high Courant number, it leads to numerical dispersion like the explicit method. On the other hand, second-order accurate semi-implicit method significantly reduces numerical dispersion. In addition to time discretization techniques, this simulator contains several space discretization methods such as first-order upstream and UMIST (University of Manchester Institute of Science and Technology) techniques. The proposed numerical simulator is suitable for easily using the different combinations of time and space discretization methods. Secondly, this study is to propose the use of semi-implicit time discretization technique with UMIST space discretization method to minimize numerical dispersion and suppress unphysical oscillation. In spite of the fact that the UMIST method suppresses unphysical oscillation, it causes a small and undesired oscillation at flood front for very large Courant number. Third objective of this study is to propose minor modification on the UMIST method to eliminate this unphysical oscillation.

Multilinear POD-DEIM model reduction for 2D and 3D semilinear systems of differential equations

<p style='text-indent:20px;'>We are interested in the numerical solution of coupled semilinear partial differential equations (PDEs) in two and three dimensions. Under certain assumptions on the domain, we take advantage of the Kronecker structure arising in standard space discretizations of the differential operators and illustrate how the resulting system of ordinary differential equations (ODEs) can be treated directly in matrix or tensor form. Moreover, in the framework of the proper orthogonal decomposition (POD) and the discrete empirical interpolation method (DEIM) we derive a two- and three-sided model order reduction strategy that is applied directly to the ODE system in matrix and tensor form respectively. We discuss how to integrate the reduced order model and, in particular, how to solve the tensor-valued linear system arising at each timestep of a semi-implicit time discretization scheme. We illustrate the efficiency of the proposed method through a comparison to existing techniques on classical benchmark problems such as the two- and three-dimensional Burgers equation.</p>

Semi-implicit Hybrid Discrete $$\left( \text {H}^T_N\right) $$ Approximation of Thermal Radiative Transfer

The thermal radiative transfer (TRT) equations form an integro-differential system that describes the propagation and collisional interactions of photons. Computing numerical solutions of the TRT equations accurately and efficiently is challenging for several reasons, the first of which is that TRT is defined on a high-dimensional phase space that includes the independent variables of time, space, and velocity. In order to reduce the dimensionality of the phase space, classical approaches such as the P\(_N\) (spherical harmonics) or the S\(_N\) (discrete ordinates) ansatz are often used in the literature. In this work, we introduce a novel approach: the hybrid discrete (H\(^T_N\)) approximation to the radiative thermal transfer equations. This approach acquires desirable properties of both P\(_N\) and S\(_N\), and indeed reduces to each of these approximations in various limits: H\(^1_N\) \(\equiv \) P\(_N\) and H\(^T_0\) \(\equiv \) S\(_T\). We prove that H\(^T_N\) results in a system of hyperbolic partial differential equations for all \(T\ge 1\) and \(N\ge 0\). Another challenge in solving the TRT system is the inherent stiffness due to the large timescale separation between propagation and collisions, especially in the diffusive (i.e., highly collisional) regime. This stiffness challenge can be partially overcome via implicit time integration, although fully implicit methods may become computationally expensive due to the strong nonlinearity and system size. On the other hand, explicit time-stepping schemes that are not also asymptotic-preserving in the highly collisional limit require resolving the mean-free path between collisions, making such schemes prohibitively expensive. In this work we develop an asymptotic-preserving numerical method that is based on a nodal discontinuous Galerkin discretization in space, coupled with a semi-implicit discretization in time. In particular, we make use of a second order explicit Runge–Kutta scheme for the streaming term and an implicit Euler scheme for the material coupling term. Furthermore, in order to solve the material energy equation implicitly after each predictor and corrector step, we linearize the temperature term using a Taylor expansion; this avoids the need for an iterative procedure, and therefore improves efficiency. In order to reduce unphysical oscillation, we apply a slope limiter after each time step. Finally, we conduct several numerical experiments to verify the accuracy, efficiency, and robustness of the H\(^T_N\) ansatz and the numerical discretizations.