A Monge-Ampère enhancement for semi-Lagrangian methods

doi:10.1016/j.compfluid.2011.01.029

Computers & Fluids

Volume 46, Issue 1, July 2011, Pages 180-185

https://doi.org/10.1016/j.compfluid.2011.01.029 Get rights and content

Abstract

Demanding the compatibility of semi-Lagrangian trajectory schemes with the fundamental Euler expansion formula leads to the Monge-Ampère (MA) nonlinear second-order partial differential equation. Given standard estimates of the departure points of flow trajectories, solving the associated MA problem provides a corrected solution satisfying a discrete Lagrangian form of the mass continuity equation to round-off error. The impact of the MA enhancement is discussed in two diverse limits of fluid dynamics applications: passive tracer advection in a steady cellular flow and in fully developed turbulence. Improvements of the overall accuracy of simulations depend on the problem and can be substantial.

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 as $x_{0} = x_{i} - \int_{t_{0}}^{t} v (x (s), s) ds,$ $ρ (x_{i}, t) = \hat{J} ρ (x_{0}, t_{0}),$ $ψ (x_{i}, t) = ψ (x_{0}, t_{0}) + \int_{T} R dt,$ where positions (x_i, t) and (x₀, t₀) are, respectively, the arrival (grid) and departure points taken along a parcel trajectory T, ρ denotes the density, $\hat{J} \equiv J^{- 1} = \det {\partial x_{0} / \partial x}$ 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).¹ In particular, for incompressible fluids considered in this paper the volume of each fluid element remains constant, upon which the relation $\hat{J} = 1$ 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 velocity $\tilde{v} \approx (Δ t)^{- 1} \int_{t_{0}}^{t} v (x (s), s) ds,$ in the discretized counterpart of (1) with the gradient of a potential ϕ $(x_{0})_{C} = x_{i} - Δ t (\tilde{v} - \nabla ϕ),$ such that the resulting set of departure points (x₀)_C satisfies $\det \{\frac{\partial (x_{0})_{C}}{\partial x}\} = 1,$ which happens to be a form of the MAE.² 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 to $A ϕ_{xx} + 2 B ϕ_{xy} + C ϕ_{yy} + E (ϕ_{xx} ϕ_{yy} - ϕ_{xy}^{2}) + D = 0,$ where subscripts x and y denote the respective partial differentiations. The coefficients of (7) depend only on partial derivatives of $\tilde{v} = (\tilde{u}, \tilde{v})$ , such that $A = 1 - Δ t {\tilde{v}}_{y}$ , $B = 1 / 2 Δ t ({\tilde{u}}_{y} + {\tilde{v}}_{x})$ , $C = 1 - Δ t {\tilde{u}}_{x}$ , $D = - {\tilde{u}}_{x} - {\tilde{v}}_{y} - Δ t ({\tilde{u}}_{y} {\tilde{v}}_{x} - {\tilde{u}}_{x} {\tilde{v}}_{y})$ 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 spheres $ψ (x, y, z) = \max (0, 1 - \sqrt{(x - x_{c})^{2} / X^{2} + (y - y_{c})^{2} / Y^{2} + (z - z_{c})^{2} / Z^{2}}),$ with X = Y = Z = 0.4π and x_c = y_c = z_c = 3π/2. The latter choice concentrates ψ near the center of the vortex at (x_c, y_c, z_c), 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)

G.I. Delzanno et al.
An optimal robust equidistribution method for two-dimensional adaptation based on Monge–Kantorovich optimization
J Comput Phys
(2008)
T. Dubois et al.
Multilevel schemes for the shallow water equations
J Comput Phys
(2005)
D.A. Knoll et al.
Jacobian-free Newton-Krylov methods: a survey of approaches and applications
J Comput Phys
(2004)
J.M. Prusa et al.
EULAG, a computational model for multiscale flows
Comput Fluids
(2008)
P.K. Smolarkiewicz et al.
A class of monotone interpolation schemes
J Comput Phys
(1992)
A. Staniforth et al.
Aspects of the dynamical core of a nonhydrostatic, deep-atmosphere, unified weather and climate-prediction model
J Comput Phys
(2008)
J. Strain
Semi-Lagrangian methods for level set equations
J Comput Phys
(1999)
D.B. Xiu et al.
A semi-Lagrangian high-order method for Navier–Stokes equations
J Comput Phys
(2001)
V. Zheligovsky et al.
The Monge-Ampère equation: various forms and numerical solution
J Comput Phys
(2010)
Y. Brenier
Polar factorization and monotone rearrangement of vector-valued functions
Comm Pure Appl Math
(1991)

R. Courant et al.

(1962)

Cited by (8)

Adaptive Isogeometric Analysis using optimal transport and their fast solvers
2024, Computer Methods in Applied Mechanics and Engineering
In this work, we devise fast solvers and adaptive mesh generation procedures based on the Monge–Ampère Equation using B-Splines Finite Elements, within the Isogeometric Analysis framework. Our approach ensures that the constructed mapping is a bijection, which is a major challenge in Isogeometric Analysis. First, we use standard B-Splines Finite Elements to solve the Monge–Ampère Equation. An analysis of this approach shows serious limitations when dealing with high variations near the boundary. In order to solve this problem, a new formulation is derived using compatible B-Splines discretization based on a discrete DeRham sequence. A new fast solver is devised in this case using the Fast Diagonalization method. Different tests are provided and show the performance of our new approach.
Optimization-based mesh correction with volume and convexity constraints
2016, Journal of Computational Physics
Citation 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.
We consider the problem of finding a mesh such that 1) it is the closest, with respect to a suitable metric, to a given source mesh having the same connectivity, and 2) the volumes of its cells match a set of prescribed positive values that are not necessarily equal to the cell volumes in the source mesh. This volume correction problem arises in important simulation contexts, such as satisfying a discrete geometric conservation law and solving transport equations by incremental remapping or similar semi-Lagrangian transport schemes. In this paper we formulate volume correction as a constrained optimization problem in which the distance to the source mesh defines an optimization objective, while the prescribed cell volumes, mesh validity and/or cell convexity specify the constraints. We solve this problem numerically using a sequential quadratic programming (SQP) method whose performance scales with the mesh size. To achieve scalable performance we develop a specialized multigrid-based preconditioner for optimality systems that arise in the application of the SQP method to the volume correction problem. Numerical examples illustrate the importance of volume correction, and showcase the accuracy, robustness and scalability of our approach.
The geometry of r-adaptive meshes generated using optimal transport methods
2015, Journal of Computational Physics
Citation 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 principles of mesh equidistribution and alignment play a fundamental role in the design of adaptive methods, and a metric tensor and mesh metric are useful theoretical tools for understanding a method's level of mesh alignment, or anisotropy. We consider a mesh redistribution method based on the Monge–Ampère equation which combines equidistribution of a given scalar density function with optimal transport. It does not involve explicit use of a metric tensor, although such a tensor must exist for the method, and an interesting question to ask is whether or not the alignment produced by the metric gives an anisotropic mesh. For model problems with a linear feature and with a radially symmetric feature, we derive the exact form of the metric, which involves expressions for its eigenvalues and eigenvectors. The eigenvectors are shown to be orthogonal and tangential to the feature, and the ratio of the eigenvalues (corresponding to the level of anisotropy) is shown to depend, both locally and globally, on the value of the density function and the amount of curvature. We thereby demonstrate how the optimal transport method produces an anisotropic mesh along a given feature while equidistributing a suitably chosen scalar density function. Numerical results are given to verify these results and to demonstrate how the analysis is useful for problems involving more complex features, including for a non-trivial time dependant nonlinear PDE which evolves narrow and curved reaction fronts.
The Monge-Ampère trajectory correction for semi-Lagrangian schemes
2014, Journal of Computational Physics
Citation 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.
Requiring that numerical estimates of the flow trajectories comply with the fundamental Euler expansion formula that governs the evolution of a volume of fluid leads to a second-order nonlinear Monge–Ampère partial differential equation (MAE). In Cossette and Smolarkiewicz (2011) [15], a numerical algorithm based on solving the MAE with an inexact Newton–Krylov solver has been developed and used to correct standard estimates of the departure points of flow trajectories provided by a classical semi-Lagrangian scheme in the context of an incompressible fluid. Here we extend the theoretical analysis of the elemental rotational and deformational motions presented in Cossette and Smolarkiewicz (2011) [15]. In particular, closed-form analytic solutions are derived for both cases that serve to illustrate the mechanics of the enhanced trajectory scheme and to address the issues of existence and uniqueness. Scalar advection shows that the MA correction improves mass conservation substantially by suppressing anomalous fluid contraction. The impact of the MA correction on complex flows is studied in the framework of an ideal magneto-fluid in which the formation of current sheets leads to topological changes and reconnection of field lines. The use of the MA correction improves the numerical stability of the solutions and prevents trajectory intersections as well as the spurious growth of the magnetic energy. 2D and 3D examples are presented and the computational performance of the solver is documented.
EULAG, a computational model for multiscale flows: An MHD extension
2013, Journal of Computational Physics
EULAG is an established high-performance computational model for simulating fluid flows across a wide range of scales and physical scenarios [Prusa et al., Comput. Fluids 37 (2008) 1193]. Historically driven by interests in simulating weather and climate processes, the numerics of EULAG are unique, owing to a synergistic blend of non-oscillatory forward-in-time MPDATA methods, robust elliptic solver, and generalized coordinate formulation enabling grid adaptivity. In this paper the numerical apparatus of an ideal magnetohydrodynamic (MHD) extension of the EULAG model is discussed, the robust workings of which have been recently revealed in global large-eddy simulations of solar magneto-convection producing solar-like magnetic cycles and dynamo action [Ghizaru et al., ApJL 715 (2010) L133; Racine et al., ApJ 735 (2011) 46]. Here, a specialized nonoscillatory forward-in-time scheme for integrating ideal anelastic MHD equations is presented in detail, and illustrated with an abstract example of magnetized three-dimensional flow in time-dependent geometry for a weak, moderate and strong magnetic field. An analysis of the model performance reveals that multiple solutions of elliptic problems do not have to imply proportionally larger computational expense.
APPLYING AN ORIENTED DIVERGENCE THEOREM TO SWEPT FACE REMAP
2023, SIAM Journal on Numerical Analysis

View all citing articles on Scopus

View full text

A Monge-Ampère enhancement for semi-Lagrangian methods

Abstract

Introduction

Section snippets

Interpretation in terms of flow shear

Passive advection in a cellular flow

Remarks

Acknowledgments

J Comput Phys

J Comput Phys

J Comput Phys

Comput Fluids

J Comput Phys

J Comput Phys

J Comput Phys

J Comput Phys

J Comput Phys

Polar factorization and monotone rearrangement of vector-valued functions

Comm Pure Appl Math