semi-Lagrangian Method Research Articles

An unconditionally stable scheme for the immersed boundary method with application in cardiac mechanics

The stability of the immersed boundary (IB) method is a challenge in simulating fluid–structure interaction problems, where time step constraints are significantly stricter than in pure fluid simulations. We propose a novel unconditionally stable scheme for the immersed boundary finite element (IBFE) method. The structure is handled implicitly and characterized by strain energy functions, rather than being modeled as fibers or membranes. Through energy estimate, we prove the unconditional stability of the fully discrete approximation in the absence of the convective term. In real simulations of cardiac mechanics problems, the time step is much larger, only limited by the Courant–Friedrichs–Lewy condition of the fluid. The novelty of this work lies in the combination of dual interpolation and distribution operators, the Jacobian-free Newton–Krylov method for solving nonlinear algebraic systems, and the semi-Lagrangian method for handling the convective term. To validate the effectiveness and accuracy of our approach, we present various benchmarks and conduct a quasi-static simulation of a three-dimensional real left ventricular model. We have shown that the numerical stability of our scheme is very robust even with much larger time step compared to conventional explicit IB methods. Our work paves the way for further works on efficient solvers of the IBFE method.

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.

Asymptotic preserving semi-Lagrangian discontinuous Galerkin methods for multiscale kinetic transport equations

We propose a new family of asymptotic preserving (AP) discontinuous Galerkin (DG) schemes for kinetic transport equations under a diffusive scaling. They are built upon a characteristics-based model reformulation proposed in Zhang et al. (2023) [27]. The reformulated model comprises a convection-diffusion-like equation for the macroscopic density and another evolution equation for the distribution function. We employ the equations in a weak form to enjoy the advantages of a DG method. We also leverage a semi-Lagrangian (SL) approach to achieve higher efficiency than a fully implicit scheme without loss of unconditional stability. A damping technique is introduced to control spurious numerical oscillations for high-order DG methods. We establish formal AP and Fourier stability analyses for fully discrete schemes in a linear scenario. Numerical experiments, covering one-dimensional to three-dimensional in space problems, confirm the accuracy, AP property, uniformly unconditional stability, and high parallel efficiency of the proposed methods.

Semi-Lagrangian approximations of the transfer operator in divergent form

Abstract The paper demonstrates two approaches to constructing monotonic difference schemes for the transfer equation in divergent form from the family of semi-Lagrangian methods: Eulerian–Lagrangian and Lagrangian–Eulerian. Within each approach, a monotonic conservative difference scheme is proposed. It is shown that within the framework of the Lagrangian–Eulerian approach, based on the use of curvilinear grids formed by the characteristics of the approximated transfer operator, it is possible to construct monotonic difference schemes of second order accuracy.

Large-Eddy Simulations with remeshed Vortex methods: An assessment and calibration of subgrid-scale models

This study assesses various subgrid-scale models within the framework of Large Eddy Simulation (LES) using a remeshed Vortex method (RVM). RVM is a semi-Lagrangian method discretizing the vorticity-velocity Navier–Stokes equations that has proven to be a stable and less dissipative alternative to more classical Eulerian methods. The subgrid-scale models are first tested on the well-known Taylor–Green Vortex case at Re=5000. Notably, the Variational Multiscale (VMS) variant of the Smagorinsky model and the Spectral Vanishing Viscosity (SVV) approaches emerge as the best-suited to the RVM, as they add diffusion to only the smallest resolved vorticity scales. Then, a stochastic uncertainty quantification analysis is conducted for both selected models, and the model coefficients are calibrated against direct numerical simulation. These coefficients are then applied to additional cases (different regimes, grid resolutions and test cases), showing the robustness of the calibration within the RVM-LES framework.

Modeling and simulation of pollution transport in the Mediterranean Sea using enriched finite element method

This paper presents a novel numerical method for simulating the transport and dispersion of pollutants in the Mediterranean sea. The governing mathematical equations consist of a barotropic ocean model with friction terms, bathymetric forces, Coriolis and wind stresses coupled to an advection–diffusion equation with anisotropic dispersion tensor and source terms. The proposed numerical solver uses a multilevel adaptive semi-Lagrangian finite element method that combines various techniques, including the modified method of characteristics, finite element discretization, coupled projection scheme based on a rotational pressure correction algorithm, and an adaptive L2-projection. The approach employs the gradient of the concentration as an error indicator for enrichment adaptations and increasing the number of quadrature points where needed without refining the mesh. The method is shown to provide accurate and efficient simulations for pollution transport in the Mediterranean sea. The proposed approach distinguishes itself from the well-established adaptive finite element methods for incompressible viscous flows by retaining the same structure and dimension of linear systems during the adaptation process.

Fourth-Order Conservative Non-splitting Semi-Lagrangian Hermite WENO Schemes for Kinetic and Fluid Simulations

AbstractWe present fourth-order conservative non-splitting semi-Lagrangian (SL) Hermite essentially non-oscillatory (HWENO) schemes for linear transport equations with applications for nonlinear problems including the Vlasov–Poisson system, the guiding center Vlasov model, and the incompressible Euler equations in the vorticity-stream function formulation. The proposed SL HWENO schemes combine a weak formulation of the characteristic Galerkin method with two newly constructed HWENO reconstruction methods. The new HWENO reconstructions are meticulously designed to strike a delicate balance between curbing numerical oscillation and introducing excessive dissipation. Mass conservation naturally holds due to the weak formulation of the semi-Lagrangian discontinuous Galerkin method and the design of the HWENO reconstructions. We apply a positivity-preserving limiter to maintain the positivity of numerical solutions when needed. Abundant benchmark tests are performed to verify the effectiveness of the proposed SL HWENO schemes.

Stabilized bases for high-order, interpolation semi-Lagrangian, element-based tracer transport

In a computational fluid model of the atmosphere, the advective transport of trace species, or tracers, can be computationally expensive. For efficiency, models often use semi-Lagrangian advection methods. High-order interpolation semi-Lagrangian (ISL) methods, in particular, can be extremely efficient, if the problem of property preservation specific to them can be addressed. Atmosphere models often use geometrically and logically nonuniform grids for efficiency and, as a result, element-based discretizations. Such grids and discretizations make stability a particular problem for ISL methods. Generally, high-order, element-based ISL methods that use the natural polynomial interpolant associated with a nodal finite-element discretization are unstable. We derive new bases having order of accuracy up to nine, with positive nodal weights, that stabilize the element-based ISL method. We use these bases to construct the linear advection operator in the property-preserving Interpolation Semi-Lagrangian Element-based Transport (Islet) method. Then we discuss key software implementation details. Finally, we show performance results for the Energy Exascale Earth System Model's atmosphere dynamical core, comparing the original and new transport methods. These simulations used up to 27,600 Graphical Processing Units (GPU) on the Oak Ridge Leadership Computing Facility's Summit supercomputer.

An error predict-correction formula of the load vector in the BSLM for solving three-dimensional Burgers’ equations

This paper aims to develop an algorithm reducing the computational cost of the backward semi-Lagrangian method for solving nonlinear convection–diffusion equations. For this goal, we introduce an error predict-correction formula (EPCF) of the load term for the Helmholtz system. The EPCF is built up to involve the same values as solving the perturbed Cauchy problem, which allows the reuse of the values. These reuses dramatically reduce the computational cost by eliminating the interpolation processes to evaluate the load term in most computational domains. The algorithm also produces a stable solution, even where the original solution has a sharp gradient, by introducing an indicator to distinguish the sharpness of the solution. Furthermore, the proposed algorithm discusses fast solvers of the discrete systems for the perturbed Cauchy problem, for which we utilize an eigenvalue decomposition theory. To show the adaptability and efficiency of the proposed algorithm, we apply it to the convection–diffusion equations such as two and three-dimensional viscous Burgers’ equations. Through several simulations, we confirm that the proposed algorithm for the load term radically reduces the computational cost, and the indicator for stiff components guarantees precision.

Semi-Lagrangian nodal discontinuous Galerkin method for the BGK model

Semi-Lagrangian nodal discontinuous Galerkin method for the BGK model

Suppressing instability in a Vlasov–Poisson system by an external electric field through constrained optimization

Maintaining the stability and shape of a plasma is a crucial task in many technological applications ranging from beam shaping to fusion energy. This is often challenging as plasma systems tend to be naturally unstable and kinetic effects can play an important role in the behavior of the instabilities. Due to the large energies involved, the primary way to control plasma systems is the application of external electromagnetic fields.In this work, we deploy a PDE-constrained optimization formulation that uses a kinetic description for plasma dynamics as the constraint. This is to optimize, over all possible controllable external electric fields, the stability of the plasma dynamics under the condition that the Vlasov–Poisson (VP) equation is satisfied. Mathematically we formulate it as a PDE-constrained optimization. For computing the functional derivative with respect to the external field in the optimization updates, the adjoint equation is derived. Furthermore, in the discrete setting, where we employ the semi-Lagrangian method as the forward solver, we also explicitly formulate the corresponding adjoint solver and the gradient as the discrete analogy to the adjoint equation and the Fréchet derivative. A distinct feature we observed of this constrained optimization is the complex landscape of the objective function and the existence of numerous local minima, largely due to the hyperbolic nature of the VP system. To overcome this issue, we utilize a gradient-accelerated genetic algorithm, leveraging the advantages of the genetic algorithm's exploration feature to cover a broader search of the solution space and the fast local convergence aided by the gradient information. We show that our algorithm obtains good electric fields that are able to maintain a prescribed profile in a beam shaping problem and uses nonlinear effects to suppress plasma instability in a two-stream configuration.

A conservative semi-Lagrangian method for inhomogeneous Boltzmann equation

The Boltzmann equation is well-known for its good predictability of behaviors of rarefied gas flows, which may not be resolved enough by continuum models. In this work, we propose a high order conservative semi-Lagrangian method for the inhomogeneous Boltzmann equation of rarefied gas dynamics. We adopt a semi-Lagrangian scheme for the convection term, which enables us to avoid CFL-type restrictions on the time step. Also, we use a fast spectral method for computation of the collision operator. In order to preserve conservative quantities, we combine a high order conservative reconstruction and a weighted optimization technique. Several numerical tests illustrate the accuracy and efficiency of the proposed method.

Parallel-in-time integration of the shallow water equations on the rotating sphere using Parareal and MGRIT

Despite the growing interest in parallel-in-time methods as an approach to accelerate numerical simulations in atmospheric modeling, improving their stability and convergence remains a substantial challenge for their application to operational models. In this work, we study the temporal parallelization of the shallow water equations on the rotating sphere combined with time-stepping schemes commonly used in atmospheric modeling due to their stability properties, namely an Eulerian implicit-explicit (IMEX) method and a semi-Lagrangian semi-implicit method (SL-SI-SETTLS). The main goal is to investigate the performance of parallel-in-time methods, namely Parareal and Multigrid Reduction in Time (MGRIT) when these well-established schemes are used on the coarse discretization levels and provide insights on how they can be improved for better performance. We begin by performing an analytical stability study of Parareal and MGRIT applied to a linearized ordinary differential equation depending on some temporal parallelization parameters, including the choice of a coarse scheme. Next, we perform numerical simulations of two standard tests in atmospheric modeling to evaluate the stability, convergence, and speedup provided by the parallel-in-time methods compared to a fine reference solution computed serially. We also conduct a detailed investigation on the influence of artificial viscosity and hyperviscosity approaches, applied on the coarse discretization levels, on the performance of the temporal parallelization. Both the analytical stability study and the numerical simulations indicate a poorer stability behavior when SL-SI-SETTLS is used on the coarse levels, compared to the IMEX scheme. With the IMEX scheme, a better trade-off between convergence, stability, and speedup compared to serial simulations can be obtained under proper parameters and artificial viscosity choices, opening the perspective of the potential competitiveness for realistic models.

Algorithms for coupled Burgers’ equations by sharing characteristic curves within BSLM

This paper introduces a new perspective of the traditional view on the velocity of each physical particle in the coupled Burgers’ equation in the backward semi-Lagrangian method (BSLM). The proposed methods reduce the number of Cauchy problems to be solved by observing a single virtual characteristic curve with a velocity. This can drastically reduce the computational cost of determining the departure point. Then, we solve the derived system reflected by the single virtual characteristic curve. Moreover, an efficient strategy for the derived linear system of equations is provided. Four examples are tested to demonstrate the adaptability and efficiency of the proposed method. The test results show that the proposed method has third- and fourth-order accuracy in time and space, respectively. In addition, compared with the existing method of solving the problem along two particles with different velocities, we confirm that the proposed method significantly reduces computational cost while maintaining accuracy well.

The dissipative Generalized Hydrodynamic equations and their numerical solution

“Generalized Hydrodynamics” (GHD) stands for a model that describes one-dimensional integrable systems in quantum physics, such as ultra-cold atoms or spin chains. Mathematically, GHD corresponds to nonlinear equations of kinetic type, where the main unknown, a statistical distribution function f(t,z,θ), lives in a phase space which is constituted by a one-dimensional position variable z, and a one-dimensional “kinetic” variable θ, actually a wave-vector, called “rapidity”. Two key features of GHD equations are first a non-local and nonlinear coupling in the advection term, and second an infinite set of conserved quantities, which prevent the system from thermalizing. To go beyond this, we consider the dissipative GHD equations, which are obtained by supplementing the right-hand side of the GHD equations with a non-local and nonlinear diffusion operator or a Boltzmann-type collision integral. In this paper, we deal with new high-order numerical methods to efficiently solve these kinetic equations. In particular, we devise novel backward semi-Lagrangian methods for solving the advective part (the so-called Vlasov equation) by using a high-order time-Taylor series expansion for the advection fields, whose successive time derivatives are obtained by a recursive procedure. This high-order temporal approximation of the advection fields is used to design new implicit/explicit Runge–Kutta semi-Lagrangian methods, which are compared to Adams–Moulton semi-Lagrangian schemes. For solving the source terms, constituted by the diffusion and collision operators, we use and compare different numerical methods of the literature.

A semi-implicit semi-Lagrangian method for simulating immersed boundary motion under high inertia and elasticity

In this paper, we present an efficient and stable fractional-step 2D immersed boundary (IB) method for solving the interaction problems between bulk fluid and elastic interface, in particular, when the fluid inertia and the interfacial elasticity are the significant factors affecting its dynamics. In myriads of real-world applications, the effects of high inertia and elasticity are dominant. So the complex fluid dynamics under such harsh conditions is an important topic in computational physics and is inherently challenging due to high computational complexity. In turn, it requires to solve the governing equations of elastic interfacial motion in an implicit manner so that more stable simulations can be performed by relaxing the Courant-Friedrichs-Lewy (CFL) condition. The contributions of the proposed approach are three folds. First, an iteration-free semi-Lagrangian method is employed in Navier-Stokes (NS) equations. Second, the elastic force acting along the interface is treated semi-implicitly in IB formulations. Both approaches improve the numerical stability associated with the high fluidic inertia and interfacial elasticity. Finally, to solve the resulting linear system, our novel idea is to transform the original 3-by-3 block matrix system into a reduced 2-by-2 block matrix system using a discrete projection operator in a staggered grid, and then explicitly represent the exact solution via the Schur complement of the Helmholtz operator. Owing to this feature, we refer to this proposed approach as reduced immersed boundary method (rIBM). We show that the two systems are equivalent in theory, whereas the conventional immersed boundary projection method (IBPM) modifies the discrete momentum equation in the original system. A series of numerical tests is conducted to confirm the stability of the rIBM using relatively larger time-step sizes, specifically with Reynolds number and inverse capillary number equal to or larger than approximately 1000. By estimating the computational time, the numerical efficiency of the proposed method is further verified in comparison with the conventional IBPM and the Crank-Nicolson scheme-based IB method. In conclusion, the proposed approach not only improves the numerical stability, but also increases the computational speed, suitable for solving more realistic problems.

Полулагранжевы аппроксимации оператора конвекции в симметричной форме

Рассмотрены два полулагранжевых численных метода для одномерного (по пространству) уравнения переноса с оператором в симметричной форме: эйлероволагранжев и лагранжево-эйлеров. Оба метода свободны от ограничения Куранта на соотношение шагов по времени и пространству. Причем во втором методе достигнут второй порядок аппроксимации для гладких решений и продемонстрировано отсутствие численной вязкости для разрывных решений. Purpose. The purpose of the study is the development and comparison of two numerical semi-Lagrangian methods with fulfillment of the conservation law at a discrete level. The approach is applied for the transport equation in the symmetric form, reflecting the law of conservation for the square of the transferred substance. The article presents the Euler – Lagrangian method, built on a rectangular difference grid that uses local values of characteristics to calculate the coefficients of difference equations. Lagrangian – Euler method is built on a spatial non-uniform grid obtained by crossing the characteristic trajectories of the equation with lines in time. Methodology. The integro-interpolation method is applied to derive approximations for the differential operator which allowed obtaining simple formulas connecting values of the grid function at the neighboring layers in time. Numerical calculations of characteristic trajectories are held by the Euler method or the Runge – Kutta method of the second order, depending on the required accuracy. Findings. Numerical methods with the mentioned properties are developed and numerically confirmed, convergence and discrete conservation laws for them are mathematically proved. The first order convergence for both time and space is proved for the Euler – Lagrange method. The second order convergence also in time and space is proved for the Lagrange – Euler method. Originality/value. The Euler – Lagrange and Lagrange – Euler methods for the numerical solution of the convection equation are developed. These methods induce differential conservation law at discrete level. The first and the second order of convergence correspondingly are mathematically proved for them. The Lagrange – Euler method has showed two improved aspects: firstly, it has greater order of convergence than the Euler – Lagrange one and secondly, it allows solving problems with the discontinuous solutions without smoothing them.

Conservative semi-Lagrangian method for continuity equation with various time steps in different parts of computational domain

The system of Navier Stokes differential equations describes gas flow near rocket nozzle. To find solution of this system in general case, scientists and engineers use numerical methods. Modern numerical methods do not allow engineers to model all features of gas flow. It happened because of sophisticated physical processes and limitations of computational hardware. So, at least where are two ways to improve it: to enlarge hardware or reduce computational complicacy. The global aim of many science investigations is to develop numerical approach with reduced computational complicacy and at the same time without loss of computational accuracy. In 1959, Aksel C. Wiin-Nielsen proposed a new and effective trajectories method for problem of numerical forecasting. In 1966, K. M. Magomedov developed similar approach (method of characteristics) for numerical modelling of space (three dimensional on space) gas flow. In 1982, O. Pironneau showed a new and effective approach for two dimensional approximation of the Navier Stokes problem. It was based on method of characteristics also. Nowadays, these methods are called semi-Lagrangian or Eulerian Lagrangian methods. They use Lagrangian nature of the transport process. To apply this advantage, scientists decompose each equation of Navier Stokes system into three parts: convective part (hyperbolic type), elliptic part and part of right-hand side. Scientists use semi-Lagrangian approach to approximate the convective parts of equations. To develop and test modern algorithms from family of semi-Lagrangian methods, we use continuity equation from the system of Navier Stokes equations. Conservative versions of semi-Lagrangian approach are based on Gauss Ostrogradsky (divergence) theorem. It allows scientists to get conservation low (balance equation) for numerical solution in norm of L1 space. Aim of our investigation is to use different time steps in different parts of computation domain. It enables us to obtain at the same time three advantages: convergence of numerical solution to exact solution, reduction of computational complicacy, implementation of conservation low (balance equation) without weight coefficients. For this purpose we decompose computational domain into two parts (subdomains) and use different time steps in them. Main complication is design algorithm in the boundaries of computational subdomains. We use one dimensional (on space) problem to demonstrate ability of developing numerical method with described advantages. Generalization of considered approach for two- or even three-dimensional cases allows engineers to model gas flow more accurately and without artificial viscosity. It is essentially important in the parts of computational domain with high level of solution gradient.

Exploring shock-capturing schemes for Particles on Demand simulation of compressible flows

In this exploratory study, we apply shock-capturing schemes within the framework of the Particles on Demand kinetic model to simulate compressible flows with mild and strong shock waves and discontinuities. The model is based on the semi-Lagrangian method, where the information propagates along the characteristics while a set of shock-capturing concepts such as the total variation diminishing and weighted essentially non-oscillatory schemes are employed to capture the discontinuities and the shock-waves. The results show that the reconstruction schemes are able to remove the oscillations at the location of the shock waves and together with the Galilean invariance nature of the Particles on Demand model, stable simulations of mild to extreme compressible benchmarks can be carried out. Moreover, the essential numerical properties of the reconstruction schemes such as their spectral analysis and order of accuracy are discussed.

Conservative stabilized Runge-Kutta methods for the Vlasov-Fokker-Planck equation

In this work, we aim at constructing numerical schemes, that are as efficient as possible in terms of cost and conservation of invariants, for the Vlasov–Fokker–Planck system coupled with Poisson or Ampère equation. Splitting methods are used where the linear terms in space are treated by spectral or semi-Lagrangian methods and the nonlinear diffusion in velocity in the collision operator is treated using a stabilized Runge–Kutta–Chebyshev (RKC) integrator, a powerful alternative of implicit schemes. The new schemes are shown to exactly preserve mass and momentum. The conservation of total energy is obtained using a suitable approximation of the electric field. An H-theorem is proved in the semi-discrete case, while the entropy decay is illustrated numerically for the fully discretized problem. Numerical experiments that include investigation of Landau damping phenomenon and bump-on-tail instability are performed to illustrate the efficiency of the new schemes.