semi-Lagrangian Algorithm Research Articles

A fully implicit, asymptotic-preserving, semi-Lagrangian algorithm for the time dependent anisotropic heat transport equation

In this paper, we extend the operator-split asymptotic-preserving, semi-Lagrangian algorithm for time dependent anisotropic heat transport equation proposed in Chacón et al. (2014) [18] to use a fully implicit time integration with backward differentiation formulas. The proposed implicit method can deal with arbitrary heat-transport anisotropy ratios χ∥/χ⊥⋙1 (with χ∥, χ⊥ the parallel and perpendicular heat diffusivities, respectively) in complicated magnetic field topologies in an accurate and efficient manner. The implicit algorithm is second-order accurate temporally and demonstrates an accurate treatment at boundary layers (e.g., island separatrices), which was not ensured by the operator-split implementation. The condition number of the resulting algebraic system is independent of the anisotropy ratio, and is inverted with preconditioned GMRES. We propose a simple preconditioner that renders the finite-dimensional linear operator compact, resulting in mesh-independent convergence rates for topologically simple magnetic fields, and convergence rates scaling as ∼(NΔt)1/4 (with N the total mesh size and Δt the timestep) in topologically complex magnetic-field configurations. We demonstrate the accuracy and performance of the approach with test problems of varying complexity, including an analytically tractable boundary-layer problem in a straight magnetic field, and a topologically complex magnetic field featuring magnetic islands with extreme anisotropy ratios (χ∥/χ⊥=1010).

An asymptotic-preserving semi-Lagrangian algorithm for the anisotropic heat transport equation with arbitrary magnetic fields

We extend the recently proposed semi-Lagrangian algorithm for the extremely anisotropic heat transport equation [Chacón et al., J. Comput. Phys., 272 (2014)] to deal with arbitrary magnetic field topologies. The original scheme (which showed remarkable numerical properties) was valid for the so-called tokamak-ordering regime, in which the magnetic field magnitude was not allowed to vary much along field lines. The proposed extension maintains the attractive features of the original scheme (including the analytical Green's function, which is critical for tractability) with minor modifications, while allowing for completely general magnetic fields. The accuracy and generality of the approach are demonstrated by numerical experiment with an analytical manufactured solution.

A performance portable implementation of the semi-Lagrangian algorithm in six dimensions

This paper describes our approach to developing a simulation software application for the fully kinetic 6D-Vlasov equation, which will be used to explore physics beyond the reduced gyrokinetic model. Simulating the fully kinetic Vlasov equation requires efficient utilization of compute and storage capabilities due to the high dimensionality of the problem. In addition, the implementation needs to be extensible regarding the physical model and flexible regarding the hardware for production runs. We start on the algorithmic background to simulate the 6-D Vlasov equation using a semi-Lagrangian algorithm. The performance portable software stack, which enables production runs on pure CPU as well as AMD or Nvidia GPU accelerated nodes, is presented. The extensibility of our implementation is guaranteed through the described software architecture of the main kernel, which achieves a memory bandwidth of almost 500 GB/s on a V100 Nvidia GPU and around 100 GB/s on an Intel Xeon Gold CPU using a single code base. We provide performance data on multiple node-level architectures discussing utilized and further available hardware capabilities. Finally, the network communication bottleneck of 6-D grid-based algorithms is quantified. A verification of physics beyond gyrokinetic theory, for the example of ion Bernstein waves, concludes the work.

Minimal Exposure Paths in Time-Varying Fields: A Semi-Lagrangian Approach

In the context of Wireless Sensor Networks (WSNs), the Minimal Exposure Path (MEP) represents an important metric to evaluate the quality of network services. Most of the current literature is concentrated on stationary sensors, while few works have addressed the existence of moving nodes. Mobile Wireless Sensor Networks (MWSNs) present a great potential for detecting invaders compared to static ones, but dealing with dynamic sensors significantly increases the coverage complexity since the overall exposure of the sensor field depends on time. Therefore, in this letter, we propose an approach to compute Minimal Exposure Paths in time-varying fields based on a control optimization method called semi-Lagrangian (SL) scheme, in such a way that an intruder will be able to penetrate the dynamic field with the lowest exposure. The SL has already been proven to reach the optimal Minimal Exposure Path (MEP) on static WSNs, but concerning dynamic nodes, the proof is much more complicated. Then, we propose a heuristics that provides convergence of the SL algorithm to a result we conjecture to be the optimal one. Results with different time-varying sensor models in obstacle-free and cluttered environments have been presented and discussed.

On the Construction of Conservative Semi-Lagrangian IMEX Advection Schemes for Multiscale Time Dependent PDEs

This article is devoted to the construction of a new class of semi-Lagrangian (SL) schemes with implicit-explicit (IMEX) Runge-Kutta (RK) time stepping for PDEs involving multiple space-time scales. The semi-Lagrangian (SL) approach fully couples the space and time discretization, thus making the use of RK strategies particularly difficult to be combined with. First, a simple scalar advection-diffusion equation is considered as a prototype PDE for the development of a high order formulation of the semi-Lagrangian IMEX algorithms. The advection part of the PDE is discretized explicitly at the aid of a SL technique, while an implicit discretization is employed for the diffusion terms. In this way, an unconditionally stable numerical scheme is obtained, that does not suffer any CFL-type stability restriction on the maximum admissible time step. Second, the SL-IMEX approach is extended to deal with hyperbolic systems with multiple scales, including balance laws, that involve shock waves and other discontinuities. A conservative scheme is then crucial to properly capture the wave propagation speed and thus to locate the discontinuity and the plateau exhibited by the solution. A novel SL technique is proposed, which is based on the integration of the governing equations over the space-time control volume which arises from the motion of each grid point. High order of accuracy is ensured by the usage of IMEX RK schemes combined with a Cauchy–Kowalevskaya procedure that provides a predictor solution within each space-time element. The one-dimensional shallow water equations (SWE) are chosen to validate the new conservative SL-IMEX schemes, where convection and pressure fluxes are treated explicitly and implicitly, respectively. The asymptotic-preserving (AP) property of the novel schemes is also studied considering a relaxation PDE system for the SWE. A large suite of convergence studies for both the non-conservative and the conservative version of the novel class of methods demonstrates that the formal order of accuracy is achieved and numerical evidences about the conservation property are shown. The AP property for the corresponding relaxation system is also investigated.

The Cauchy-Lagrange method for 3D-axisymmetric wall-bounded and potentially singular incompressible Euler flows

We treat the incompressible, and axisymmetric Euler equations for a three-dimensional cylindrical domain with boundaries. The equations are solved by the novel Cauchy-Lagrange algorithm (CLA), which uses the time-analyticity of the Lagrangian trajectories of an incompressible Euler flow and computes the time-Taylor coefficients of the Lagrangian map via recursion relations. This semi-Lagrangian algorithm uses a pseudo-spectral type approach in space by approximating the flow fields by Chebyshev-Fourier polynomials. New methods are presented to solve the resulting Poisson problems directly for their second-order space derivatives. The flow fields, known on the Lagrangian trajectories after one time-step, are interpolated back onto the Eulerian grid to start a new recursion cycle. The time-step is only limited by the radius of convergence and, thus, independent of any Courant-Friedrichs-Lewy (CFL) condition. This allows to advance the flow with larger time-steps independently of the mesh. Stationary and swirl-free flows are used to thoroughly test our implementation for the given geometry. In this work, our ultimate goal is to apply the CLA to a flow, which might develop finite time singularities, resulting in a loss of smoothness. This demands either 3D, or 2D with swirl, and may well require the presence of solid boundaries, as indicated in recent numerical work (Luo and Hou (2014) [45]).

A high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangian drift–diffusion models coupled with Eulerian discontinuous spectral element method

The explicit semi-Lagrangian method for solution of Lagrangian transport equations as developed in Natarajan and Jacobs (2020) is adopted for the solution of stochastic differential equations that is consistent with Discontinuous Spectral Element Method (DSEM) approximations of Eulerian conservation laws. The method extends the favorable properties of DSEM that include its high-order accuracy, its local and boundary fitted properties and its high performance on parallel platforms for the concurrent Monte-Carlo, semi-Lagrangian and Eulerian solution of a class of time-dependent problems that can be described by coupled Eulerian–Lagrangian formulations. Such formulations include the probabilistic models used for the simulation of chemically reacting turbulent flows or particle-laden flows. Consistent with an explicit, DSEM discretization, the semi-Lagrangian method seeds particles at Gauss quadrature collocation nodes within a spectral element. The particles are integrated explicitly in time according to a drift velocity and a Wiener increment forcing and form the nodal basis for an advected interpolant. This interpolant is mapped back in a semi-Lagrangian fashion to the Gauss quadrature points through a least squares fit using constraints for element boundary values. Stochastic Monte-Carlo samples are averaged element-wise on the quadrature nodes. The stable explicit time step Wiener increment is sufficiently small to prevent particles from leaving the element’s bounds. The semi-Lagrangian method is hence local and parallel and does not have the grid complexity, and parallelization challenges of the commonly used Lagrangian particle solvers in particle-mesh methods for solution of Eulerian–Lagrangian formulations. Formal proof is presented that the semi-Lagrangian algorithm evolves the solution according to the Eulerian Fokker–Planck equation. Numerical tests in one and two dimensions for drift–diffusion problems show that the method converges exponentially for constant and non-constant advection and diffusion velocities.

The moving mesh semi-Lagrangian MMSISL method

We introduce a novel location-based moving mesh algorithm MMSISL in which the arrival points in the Semi-Implicit Semi-Lagrangian (SISL) algorithm are located by using an equidistribution strategy. This algorithm gives a natural coupling between moving mesh methods and SISL methods. It involves little extra cost in implementation as it exploits the interpolation methods already embedded in the SISL algorithm. We apply this method to a number of partial differential equation problems in one-dimension, each of which have sharply defined features. We show that using MMSISL leads to a markedly improved performance over fixed mesh methods, with significantly reduced errors. We also show that unlike many adaptive schemes, no issues arise in the MMSISL algorithm from a CFL condition imposed restriction on the time step.

A semi-Lagrangian algorithm in policy space for hybrid optimal control problems

The mathematical framework of hybrid system is a recent and general tool to treat control systems involving control action of heterogeneous nature. In this paper, we construct and test a semi-Lagrangian numerical scheme for solving the Dynamic Programming equation of an infinite horizon optimal control problem for hybrid systems. In order to speed up convergence, we also propose and analyze an acceleration technique based on policy iteration. Finally, we validate the approach via some numerical tests in low dimension.

Numerical solution of the Navier-Stokes equations using the Path Tubes method

The present work addresses an extension of the Path Tubes method for solving the time-dependent Navier-Stokes equations for an incompressible Newtonian fluid. The approach used is a physically intuitive methodology whose formulation is based on the theoretical foundations of the mechanics of continuous media. This version of the Path Tubes method draws on a semi-Lagrangian time discretization that employs the Reynolds’s transport theorem and a localization strategy. This time discretization can be seen as a transformation that acts on the Navier-Stokes equations, transforming this classical nonlinear model into linear partial differential equations of the (essentially) parabolic type.The result is an implicit semi-Lagrangian algorithm that allows the use of classical schemes for spatial discretization such as central-difference formulas, without the need to use upwind techniques or high-order corrections for time derivatives. The Path Tubes method was implemented through parallel computing. For this, we use a computer equipped with shared-memory multiprocessors and the OpenMP software. After intensive numerical tests and using different values of Reynolds numbers typical of advection-dominated flows, the proposed method proved to be accurate and able to work with coarse grids.

Application of the Path Tubes Method to the Navier-Stokes Equations

This work deals with an extension of the Path Tubes method for the solution of the time-dependent Navier-Stokes equations for an incompressible Newtonian fluid. Departing from a physically intuitive methodology based on the theoretical basis of the mechanics of continuous media, a robust numerical technique is obtained. This version of the Path Tubes method draws on a semi-Lagrangian time-discretization that employs the Reynolds’ transport theorem, and a localization approach, to establish an implicit semi-Lagrangian algorithm that allows the use of classical schemes for spatial discretization, such as central-difference formulas, without the need to use upwind techniques, or high-order corrections for time derivatives. Some of the extensive numerical tests are shown herein, in particular for Reynolds’ numbers typical of advection dominated flows. The tests show the method is accurate, even for coarse grids.

Toward a PV-Based Algorithm for the Dynamical Core of Hydrostatic Global Models

Abstract The diabatic contour-advective semi-Lagrangian (DCASL) algorithms previously constructed for the shallow-water and multilayer Boussinesq primitive equations are extended to multilayer non-Boussinesq equations on the sphere using a hybrid terrain-following–isentropic (σ–θ) vertical coordinate. It is shown that the DCASL algorithms face challenges beyond more conventional algorithms in that various types of damping, filtering, and regularization are required for computational stability, and the nonlinearity of the hydrostatic equation in the σ–θ coordinate causes convergence problems with setting up a semi-implicit time-stepping scheme to reduce computational cost. The prognostic variables are an approximation to the Rossby–Ertel potential vorticity Q, a scaled pressure thickness, the horizontal divergence, and the surface potential temperature. Results from the DCASL algorithm in two formulations of the σ–θ coordinate, differing only in the rate at which the vertical coordinate tends to θ with increasing height, are assessed using the baroclinic instability test case introduced by Jablonowski and Williamson in 2006. The assessment is based on comparisons with available reference solutions as well as results from two other algorithms derived from the DCASL algorithm: one with a semi-Lagrangian solution for Q and another with an Eulerian grid-based solution procedure with relative vorticity replacing Q as the prognostic variable. It is shown that at intermediate resolutions, results comparable to the reference solutions can be obtained.

3D conservative cascade semi-Lagrangian transport scheme using reduced latitude–longitude grid (CCS-RG)

The article presents conservative cascade scheme using reduced grid (CCS-RG) - a finite volume semi-Lagrangian (FVSL) algorithm for numerical solution of 3D advection equation on the sphere. The CCS-RG is based on the cascade approach that allows to split 3D advected quantity mass remapping between the Lagrangian and Eulerian grids onto 3 1D remappings. Another feature is using a quasi-uniform spaced reduced latitude-longitude grid on the sphere. The scheme is inherently mass-conservative, monotonic and stable with large time-steps. The CCS-RG performance is tested with the advection cases from the Dynamical Core Model Intercomparison Project test suite. The impact of the reduced grids and two different monotonic filters on the solution is investigated. The comparison with the results obtained by other schemes proves competitive accuracy of CCS-RG. In most cases, the numerical solution obtained with various reduced grids is as accurate as the regular grid solution.

Linear perspective shape-from-shading method with two images

Shape-from-shading (SFS) is one of the important approaches of 3-D surface reconstruction in computer vision. Since reflectance map equation in SFS is a nonlinear partial differential equation (PDE) with two unknown variables, SFS with one image is ill-posed in mathematical sense. A linear perspective SFS method with two images is proposed to deal with the problem. We assume that two images with different light source directions are captured firstly. Orthogonal projection is not as accurate as perspective one to simulate imaging processes. Two reflectance map equations are established based on the Lambertian model under perspective projection, and the equations are further transformed into one linear PDE. Then the iterative semi-Lagrangian algorithm is used to approximate the solution. Finally, 3-D height values of pixel points in imaging planes are solved by the numerical interpolation method. Experimental results of both hemisphere and complex surfaces show that the proposed method can reconstruct surfaces accurately.

An asymptotic-preserving semi-Lagrangian algorithm for the time-dependent anisotropic heat transport equation

We propose a semi-Lagrangian numerical algorithm for a time-dependent, anisotropic temperature transport equation in magnetized plasmas in regimes with negligible variation of the magnitude of the magnetic field B along field lines. The approach is based on a formal integral solution of the parallel (i.e., along the magnetic field) transport equation with sources. While this study focuses on a Braginskii (local) heat flux closure, the approach is able to accommodate nonlocal parallel heat flux closures as well. The numerical implementation is based on an operator-split formulation, with two straightforward steps: a perpendicular transport step (including sources), and a Lagrangian (field-line integral) parallel transport step. Algorithmically, the first step is amenable to the use of modern iterative methods, while the second step has a fixed cost per degree of freedom (and is therefore algorithmically scalable). Accuracy-wise, the approach is free from the numerical pollution introduced by the discrete parallel transport term when the perpendicular to parallel transport coefficient ratio χ⊥/χ∥ becomes arbitrarily small, and is shown to capture the correct limiting solution when ϵ=χ⊥L∥2χ∥L⊥2→0 (with L∥, L⊥ the parallel and perpendicular diffusion length scales, respectively). Therefore, the approach is asymptotic-preserving. We demonstrate the performance of the scheme with several numerical experiments with varying magnetic field complexity in two dimensions, including the case of heat transport across a magnetic island in cylindrical geometry in the presence of a large guide field.

Accurate Semi-Lagrangian Time Stepping for Stochastic Optimal Control Problems with Application to the Valuation of Natural Gas Storage

We describe a semi-Lagrangian time-stepping algorithm for a particular class of stochastic optimal control problems, applicable to storage valuation problems. The discretization in time uses a semi-Lagrangian approach based on Strang splitting, and convergence to the unique viscosity solution is established by appealing to the framework of Barles and Souganidis [Asymptotic Anal., 4 (1991), pp. 271--283]. The approach is illustrated in the context of a natural gas storage setting. A fully discrete approximation for the storage valuation problem using a Fourier-cosine method is described, and second-order convergence demonstrated, for pure-diffusion and jump-diffusion models.

On contour crossings in contour-advective simulations – part 2 – analysis of crossing errors and methods for their prevention

This is the second of two papers devoted to the analysis of contour crossing errors that occur in contour-advective simulations of fluid motion, where either vorticity or potential vorticity is represented by contours. We begin with a detailed discussion on some of the potential mechanisms for contour crossing. Past work has suggested that the formation of contour crossings is due to inadequate spatial resolution of contours [1]. The implementation of two schemes for preventing contour crossings within the framework of the Contour-Advective Semi-Lagrangian (CASL) algorithm is detailed here. We then present an analysis of contour crossing errors in simulations of quasigeostrophic turbulence on the f-plane and the quasigeostrophic motion of an initially circular vortex patch on the β-plane using the algorithm detailed in Part 1. We find that in general individual crossings occur at scales smaller than the inversion grid scale on which velocity is calculated, but at scales larger than that of the surgical scale that defines the smallest resolved features (vorticity) of a flow. If the resolution of a quasigeostrophic turbulence simulation on the f-plane is increased by doubling the number of grid points in each coordinate direction used in the calculation of the velocity field, then the total area in error due to contour crossings remains unchanged; a smaller number of crossings introducing larger scale area errors is replaced by a greater number of smaller local errors. Uniformly increasing the density of nodes along all contours and placement of nodes at points of close approach on contours are both effective methods for limiting contour crossings.

Revisiting Vacillations in Shallow-Water Models of the Stratosphere Using Potential-Vorticity-Based Numerical Algorithms

Abstract Polar vortex vacillations are investigated using long-term simulations of potential-vorticity (PV)-based shallow-water (SW) models for the stratosphere. In the models examined, mechanical forcing is applied through a time-independent topography mimicking tropospheric excitation of the stratosphere. Thermal forcing is applied through a linear relaxation of the mass field to a time-independent equilibrium state mimicking the radiative relaxation taking place in the stratosphere. The SW equations in the PV, velocity divergence, and acceleration divergence representation are solved for a range of resolutions using the “diabatic contour-advective semi-Lagrangian” (DCASL) algorithm and a standard pure semi-Lagrangian (SL) algorithm. Using very different numerical algorithms enables the determination of the degree of numerical sensitivity and the properties of the vacillations with much greater accuracy than in previous related studies. The focus here is on the Lagrangian or material evolution of the polar vortex. The authors examine quasi-Lagrangian diagnostics based on equivalent latitude, the mass enclosed by PV contours, and the terms involved in its time evolution. The PV field forms the basis for calculating quasi-Lagrangian diagnostics. Variations in the mass enclosed by a PV contour are associated with nonconservative processes such as diabatic heating, friction, and irreversible small-scale mixing. Generally, the mass of the polar vortex increases under the action of diabatic mass fluxes, whereas it decreases under the action of dissipative mass fluxes. The results herein differ from previous results reported at T42 resolution by Rong and Waugh in which a spectral transform algorithm is used to solve the SW equations in a vorticity–divergence–mass representation, and in which dissipation is provided by explicitly damping vorticity using hyperdiffusion. Except for the first large-amplitude oscillation, there is little sign of a clear, systematic phase shift between the dissipative and diabatic mass fluxes across the edge of the polar vortex, as proposed by Rong and Waugh as the main mechanism responsible for the vacillations. Concomitant with the absence of a phase shift, the vacillations tend to decay and occur intermittently. Rather than a phase shift, inherent fluctuations in both the diabatic and mass fluxes across the edge of the polar vortex appear to be responsible for the vacillations.

Developing new approaches for the path tubes method

The path tubes (PT) method for advection equations is a semi-Lagrangian algorithm recently introduced in the literature. PT method is simple, physically intuitive, backward-in-time and mass-conserving algorithm, whose formulation is based on the fundaments of the mechanics of continuous media, in the light of Lagrangian methodology. In present work, we describe new approaches for PT methodology, including an absolutely accurate semi-Lagrangian scheme for 1D problems and a new and attractive discretization for multidimensional problems. In addition, we show that PT method is unconditionally stable.

A massively parallel semi-Lagrangian algorithm for solving the transport equation

The scalar transport equation underpins many models employed in science, engineering, technology and business. Application areas include, but are not restricted to, pollution transport, weather forecasting, video analysis and encoding (the optical flow equation), options and stock pricing (the Black-Scholes equation) and spatially explicit ecological models. Unfortunately finding numerical solutions to this equation which are fast and accurate is not trivial. Moreover, finding such numerical algorithms that can be implemented on high performance computer architectures efficiently is challenging. In this paper the authors describe a massively parallel algorithm for solving the advection portion of the transport equation. We present an approach here which is different to that used in most transport models and which we have tried and tested for various scenarios. The approach employs an intelligent domain decomposition based on the vector field of the system equations and thus automatically partitions the computational domain into algorithmically autonomous regions. The solution of a classic pure advection transport problem is shown to be conservative, monotonic and highly accurate at large time steps. Additionally we demonstrate that the algorithm is highly efficient for high performance computer architectures and thus offers a route towards massively parallel application.