A Monge-Ampère enhancement for semi-Lagrangian methods
Introduction
The research on the Monge-Ampère equation (hereafter the MA equation or MAE) extends over more than two centuries [13]. Nowadays, MAE is at the core of many theoretical and computational applications deriving from functional and geometric analysis, geophysical fluid dynamics, cosmology, optimization and imaging technology; see [1], [12], [3], [4], [22], [9] and references therein. Here, we arrive at MAE while examining the impact of enforcing the compatibility of the trajectory schemes with mass/volume continuity in semi-Lagrangian (SL) integrations of equations governing fluid dynamics [8], [15], [20], [21]. These equations can be written compactly aswhere positions (xi, t) and (x0, t0) are, respectively, the arrival (grid) and departure points taken along a parcel trajectory T, ρ denotes the density, is the inverse flow Jacobian and ψ symbolizes a specific fluid variable (e.g., temperature or a velocity component) with R indicating its associated forcing.
Although SL schemes are established in CFD, and used operationally in numerical weather prediction [18], approximations to the trajectory integrals (1) are typically incompatible with the fundamental Euler expansion formula – dlnJ/dt = ∇ · v, with d/dt denoting the total derivative along T – that underlies the Lagrangian form of the mass continuity Eq. (2).1 In particular, for incompressible fluids considered in this paper the volume of each fluid element remains constant, upon which the relation becomes a necessary condition for the scheme compatibility.
On the basis that vortical motions do not induce any change in the volume of fluid elements, we correct the estimated path-mean velocityin the discretized counterpart of (1) with the gradient of a potential ϕsuch that the resulting set of departure points (x0)C satisfieswhich happens to be a form of the MAE.2 An exact-projection-type solver capable of delivering round-off error accurate solutions to (6) has been developed, and its technical exposition will be given elsewhere. Here it is used to show the impact of (5) in two diverse applications on periodic domains: a 3D passive tracer advection in a steady cellular flow, and a freely evolving 2D turbulence – archetypes of pollution transport and large-scale atmospheric circulations, respectively. In Section 2 we discuss the derivation of MAE, its physical interpretation in terms of flow rotation and deformation, and the implementation of the MAE solver in the context of a semi-Lagrangian scheme. Numerical results are presented in Section 3 and remarks in Section 4 conclude the paper.
Section snippets
Interpretation in terms of flow shear
Here we present the 2D case of (6), the archetype of which is considered in Section 4.3 of [2]. For any Δt > 0, the insertion of (5) into (6) leads towhere subscripts x and y denote the respective partial differentiations. The coefficients of (7) depend only on partial derivatives of , such that , , , and E = Δt.
Passive advection in a cellular flow
Fig. 2 shows isosurfaces of tracer distributions ψ(x, y, z) = c, identifiable with the boundaries of fluid elements, after they have been advected with the flow deriving from the stream-function Ψ(x, y, z) = (0, 0, sin(x) sin(y)); (x, y, z) ∈ [0, 2π] × [0, 2π] × [0, 2π]. The initial distribution is given by the family of concentric sphereswith X = Y = Z = 0.4π and xc = yc = zc = 3π/2. The latter choice concentrates ψ near the center of the vortex at (xc, yc, zc), where flow
Remarks
Together, the analysis in terms of the elemental flows carried in Section 2 and the cellular flow advection experiments support the importance of enforcing the compatibility of semi-Lagrangian trajectory schemes with the Euler expansion formula. The advection experiments have shown that both the shape of the fluid elements and the total mass are better preserved by the MA-enhanced algorithm. Furthermore, the experiments with freely evolving 2D turbulence have documented the better
Acknowledgments
The discussions with Paul Charbonneau are gratefully acknowledged. Comments from two anonymous referees helped to improve the presentation. This work was supported in part by the DOE award DE-FG02-08ER64535. The National Center for Atmospheric Research is sponsored by the National Science Foundation.
References (22)
- et al.
An optimal robust equidistribution method for two-dimensional adaptation based on Monge–Kantorovich optimization
J Comput Phys
(2008) - et al.
Multilevel schemes for the shallow water equations
J Comput Phys
(2005) - et al.
Jacobian-free Newton-Krylov methods: a survey of approaches and applications
J Comput Phys
(2004) - et al.
EULAG, a computational model for multiscale flows
Comput Fluids
(2008) - et al.
A class of monotone interpolation schemes
J Comput Phys
(1992) - et al.
Aspects of the dynamical core of a nonhydrostatic, deep-atmosphere, unified weather and climate-prediction model
J Comput Phys
(2008) Semi-Lagrangian methods for level set equations
J Comput Phys
(1999)- et al.
A semi-Lagrangian high-order method for Navier–Stokes equations
J Comput Phys
(2001) - et al.
The Monge-Ampère equation: various forms and numerical solution
J Comput Phys
(2010) Polar factorization and monotone rearrangement of vector-valued functions
Comm Pure Appl Math
(1991)
Cited by (8)
Adaptive Isogeometric Analysis using optimal transport and their fast solvers
2024, Computer Methods in Applied Mechanics and EngineeringOptimization-based mesh correction with volume and convexity constraints
2016, Journal of Computational PhysicsCitation Excerpt :First, provided the problem has a solution, the computed NLP solution provably satisfies the volume constraint (3) in exact arithmetic, and to near machine precision in practice. This is an important advantage over the methods from the first category and the MA trajectory correction [17]. Second, the NLP solution is locally optimal; if the original deformed mesh is chosen as the initial guess for the NLP solver (in addition to being part of the objective function definition), then the NLP solution is a constraint-consistent mesh that is in some sense closest to the original constraint-inconsistent mesh.
The geometry of r-adaptive meshes generated using optimal transport methods
2015, Journal of Computational PhysicsCitation Excerpt :This requires solving an associated scalar Monge–Ampère equation and constructing the mesh mapping from the gradient of its solution. These methods have the potential advantages of being robust, flexible, and cheap to implement, for both two and three dimensional problems, particularly CFD type problems [15,16]. They also have certain very desirable properties, such as an absence of mesh tangling and an inheritance of self-similar behaviour in the solution [17,10].
The Monge-Ampère trajectory correction for semi-Lagrangian schemes
2014, Journal of Computational PhysicsCitation Excerpt :Such properties are, for instance, conservation laws, monotonicity, positive-definiteness of the scheme, total-variation diminishing constraints, or some form of compatibility [54]. In the brief proceedings paper [15], the Monge–Ampère (hereafter, MA) trajectory correction for semi-Lagrangian (SL) integration schemes was introduced as a means to enforce a type of compatibility based on the Lagrangian form of the mass continuity equation in the context of incompressible fluids. The execution of the MA correction requires solving a diagnostic nonlinear second-order PDE, the Monge–Ampère equation (MAE), to some approximation at each model time step.
EULAG, a computational model for multiscale flows: An MHD extension
2013, Journal of Computational PhysicsAPPLYING AN ORIENTED DIVERGENCE THEOREM TO SWEPT FACE REMAP
2023, SIAM Journal on Numerical Analysis