An optimal robust equidistribution method for two-dimensional grid adaptation based on Monge–Kantorovich optimization

https://doi.org/10.1016/j.jcp.2008.07.020Get rights and content

Abstract

A new cell-area equidistribution method for two-dimensional grid adaptation, based on Monge–Kantorovich optimization (or Monge–Kantorovich optimal transport), is presented. The method is based on a rigorous variational principle, in which the L2 norm of the grid displacement is minimized, constrained locally to produce a prescribed positive-definite cell volume distribution. The procedure involves solving the Monge–Ampère equation: A single, nonlinear, elliptic scalar equation with no free parameters, and with proved existence and uniqueness theorems. We show that, for sufficiently small grid displacement, this method also minimizes the mean grid-cell distortion, measured by the trace of the metric tensor. We solve the Monge–Ampère equation numerically with a Jacobian-Free Newton–Krylov method. The ellipticity property of the Monge–Ampère equation allows multigrid preconditioning techniques to be used effectively, delivering a scalable algorithm under grid refinement. Several challenging test cases demonstrate that this method produces optimal grids in which the constraint is satisfied numerically to truncation error. We also compare this method to the well known deformation method [G. Liao, D. Anderson, Appl. Anal. 44 (1992) 285]. We show that the new method achieves the desired equidistributed grid using comparable computational time, but with considerably better grid quality than the deformation method.

Introduction

Equidistribution has traditionally been a fundamental guiding principle in grid generation, as is evidenced by the ample literature on the subject (see e.g. [1], [2], [3], [4] and references therein). The problem is very simply posed: Generate a grid that equidistributes a given quantity along an arc (1D), a surface (2D), or a volume (3D). (Hereafter, we will refer to equidistribution with respect to volumes for arbitrary dimensionality.) The concept is most attractive due to its conceptual simplicity. Furthermore, in the context of error equidistribution, a rigorous connection exists between error equidistribution and minimization of total error [5], [6].

In one dimension, equidistribution determines the grid uniquely [1], [2], [7], [8]. This follows because, in 1D, only one unknown exists per cell, and it can be determined uniquely by specifying the Jacobian of the transformation at each cell. However, in two or more dimensions, equidistribution itself is not sufficient to determine the grid uniquely. There are many possible grids that satisfy a given equidistribution principle, and therefore a possibility arises to select an optimal grid in some reasonable sense.

Historically, many approaches have been developed to determine good quality grids in two and higher dimensions. Variational methods have received a great deal of attention, as they provide a solid mathematical foundation, leading to Euler–Lagrange equations that govern the generation of the grid. There are several very good reviews on the subject [1], [2], [7], [8]. Pioneering work by Winslow [9] used the “smoothness” integral (a measure of the trace of the metric tensor) to determine a set of Laplace equations that govern the generation of the grid (equipotential method). This method was very successful in allowing the modeler to generate grids for complicated boundaries. In a later development [10], Winslow proposed the variable diffusion method, in order to have control of the properties of the grid within the domain. Brackbill and Saltzman [11] used a cost function consisting of a combination of smoothness, orthogonality, and volume variation integrals to define a grid generation equation that incorporates properties of each. In this formulation, a term related to equidistribution is introduced by the volume variation integral. However, local equidistribution was not achieved due to competition among the different cost functions. Dvinsky [12] pioneered the use of harmonic maps for grid generation, and they have since been explored by a number of authors [4], [13], [14], [15]. Harmonic maps can also be derived from variational principles, and they are very attractive because, in 2D under certain conditions, the existence of a solution is guaranteed [16], [17]. However, they have important drawbacks, as we outline in the next paragraph. More recently, Huang [18] revisited the issue of variational equidistribution by providing rigorous integral measures of isotropy (or smoothness) and uniformity (or equidistribution). Again, equidistribution appears as a fundamental ingredient of the grid generation strategy.

Despite the attractiveness of these variational approaches – their mathematical soundness and the reasonable quality of the resulting grids – they suffer from various drawbacks. In particular, some require user-provided parameters to decide the relative contributions of the globally averaged terms in the cost function, and, when these parameters are chosen poorly, this may lead to mathematically ill-posed problems [19]. In addition, they result in as many Euler–Lagrange equations as dimensions considered, and these are strongly coupled and very nonlinear. They are therefore difficult to solve numerically (although, in the context of harmonic functionals, there has been recent progress employing state-of-the-art nonlinear algorithms [20]). The major disadvantage of some of these approaches is the following: Because the global grid property integrals compete against each other, the grid never truly satisfies any constraint, including equidistribution, to any predictable accuracy.

To resolve some of these issues, elliptic grid generation methods were first proposed by Thompson et al. [21], which evolved from the earlier work by Winslow [9]. In these approaches, non-homogeneous terms are added to Winslow’s equipotential method. When properly specified, such terms allow good control of the properties of the grid. In [22], the source terms are found from a least squares (variational) fit of the inverse Jacobian matrix of the transformation to a target matrix with the desired properties. In Refs. [3], [23], evidence that an elliptic approach can be made equivalent to an equidistribution principle by a proper choice of the non-homogeneous terms is presented. The attractiveness of these methods is that the equations remain essentially elliptic, and are therefore fairly tractable algorithmically. However, they still require as many coupled, nonlinear elliptic equations as the dimensionality of the problem, there is no rigorous existence and uniqueness theory, and for the most part they lack the mathematical soundness of variational principles.

Several authors have attempted to generalize the concept of grid equidistribution directly to multiple dimensions. One approach is to consider equidistribution along one-dimensional arcs in the multi-dimensional domain [3], [24], [25]. This method has the advantage that the task of multidimensional grid generation may be decomposed in a series of 1D equidistribution steps along coordinate arcs [3], [26], [27]. However, it has been shown [3], [24] that the concept of arc equidistribution can only be satisfied locally in the domain, not globally, and that it generates fairly poor-quality grids (and may even fold the grid [24]).

The case for the need of cell-volume equidistribution to fix such smoothness problems has been made by various authors [28], [23]. As we have argued earlier, in dimensions greater than unity equidistribution does not guarantee a unique solution. This implies that there is room for grid optimization. Recently, Kania [23] derived a set of non-homogeneous terms for Thompson’s method that achieves volume equidistribution, and demonstrated that the generated grids are of good quality (although no a posteriori measure is given of how accurate the equidistribution principle is satisfied by the generated grids). Liao and Anderson [29] proposed an ODE-based equidistribution approach based on the work of Moser [30], [31]. By suitably defining a flow velocity and an accompanying set of ODEs, it was demonstrated that the approach leads to an equidistributed grid. Then, based on ODE theory, they demonstrate that the solution obtained by the procedure exists and is unique, once the flow is prescribed. However, there is great latitude in choosing the flow, and for a fixed flow there is no evidence that the resulting grid is optimal in any sense.

In this paper, we propose a new approach for cell-volume equidistribution, based on Monge–Kantorovich optimization [32], [33]. This method is based on a constrained minimization approach. Instead of minimizing a quantity consisting of a grid quality measure plus an equidistribution measure, this method involves minimizing a grid quality measure constrained locally by the equidistribution principle. This constraint is enforced by a local Lagrange multiplier. In this fashion, the method chooses the optimal grid which is compatible with the equidistribution principle. The minimization procedure results in a single, nonlinear, elliptic equation for the Lagrange multiplier with no tunable parameters, the Monge–Ampère equation [34]. This equation has been shown (see e.g. [35]) to have a unique solution in 2D and 3D. Our variational method achieves equidistribution up to truncation error. This is unlike existing elliptic or variational approaches which, besides requiring as many equations as the dimensionality of the problem, in general do not enforce equidistribution locally. Furthermore, in this work, we also establish a connection between the Monge–Kantorovich approach and the minimization of grid cell distortion (smoothness), which is a very desirable property in grid generation (see e.g. [11], [24]).

Our approach combines the advantages of both variational and elliptic grid generation approaches. The variational character of our method is obvious, as it stems from a minimization procedure. The leading term of the resulting Euler–Lagrange equation is a Laplace operator (Section 2.2), and the linearized full operator is elliptic in nature (Section 3.2). Therefore, the method is suitable for modern, fast nonlinear solvers for elliptic equations. Specifically, we will demonstrate in this study the effectiveness of multigrid-preconditioned Newton–Krylov methods.

We would like to point out that the Monge–Ampère equation has been discussed before in the context of grid generation in Ref. [36], where its suitability for blow-up problems (with a developing point singularity) was shown. However, in that reference, the Monge–Ampère equation is parabolized to obtain an approximate solution (vs. the scalable, fully nonlinear algorithm proposed here). Further, there was no discussion of grid optimality in any sense.

The remainder of this paper is organized as follows. In Section 2 we formulate the problem and introduce our equidistribution approach based on Monge–Kantorovich optimization, and the resulting Monge–Ampère equation. In Section 3 we discuss certain properties of the Monge–Kantorovich approach, such as ellipticity and the connection with the minimization of grid distortion (maximization of grid smoothness). The numerical implementation of the equidistribution PDE is briefly discussed in Section 4. In Section 5, the Monge–Kantorovich approach is tested with several challenging examples. A comparison is also made with the deformation method [29] and with a method that minimizes grid-cell distortion. All the tests demonstrate the effectiveness and robustness of the Monge–Kantorovich approach in achieving optimally equidistributed grids. Conclusions are drawn in Section 6. The Appendix briefly reviews the deformation method [29].

Section snippets

Prescribing the Jacobian

The problem involves finding a one-to-one transformation in physical space according to a prescribed transformation Jacobian (or density or monitor function), i.e. to generate an adaptive grid with prescribed volumes. In this study, we focus on the 2D case. An example of the applicability of the method to 3D is presented in Ref. [37]. Let XR2 be a bounded domain with boundary X. We define a two-dimensional coordinate transformation in physical space between the coordinates of an initial grid x

Properties of Monge–Kantorovich optimization

We now proceed to establish certain important properties of the Monge–Kantorovich approach and the Monge–Ampère equation, Eq. (21). These include the relation with grid distortion (grid smoothness) and ellipticity.

Numerical implementation

In this section, we discuss the details of the numerical implementation of the Monge–Kantorovich approach, Eq. (21). Here and in Section 5, for simplicity, we will consider X=Ξ and the identity mapping between the logical space Ξ and the initial grid in physical space: x=ξ. This implies J=1.

Results

In this Section, we will apply the Monge–Kantorovich approach to several challenging tests. First, we will solve the inverse problem, Eq. (22), to compare the Monge–Kantorovich approach to the optimal distortion method (inverse approach), in order to test the connection between the two approaches established in Section 3.1 for an example where εO(1).

Second, we will compare the performance of the Monge–Kantorovich approach against the deformation method of Liao and Anderson [29]. (Again we will

Conclusions

We have presented an effective and robust approach for cell-area equidistribution in 2D. The new method, which is based on Monge–Kantorovich optimization, is obtained by a minimization of the L2 norm of the grid displacement, constrained to satisfy a given distribution of cell volumes.

This method is based on a single, nonlinear, scalar equation, the Monge–Ampère equation. This equation has no free parameters, and solutions are known to exist and be unique. We have solved the Monge–Ampère

Acknowledgment

The authors are very grateful to Jerry Brackbill, Darryl Holm, Alex Dragt, Brent Wohlberg, and Richard Chartrand for stimulating discussions. GLD wishes to thank Caterina Baravalle for constant support. This work was funded by the Laboratory Directed Research and Development program (LDRD), US Department of Energy Office of Science, Office of Fusion Energy Sciences, under the auspices of the National Nuclear Security Administration of the US Department of Energy by Los Alamos National

References (61)

  • J.F. Thompson et al.

    Automatic numerical generation of body fitted curvilinear coordinate system for field containing any number of arbitrary two dimensional bodies (appl. to airfoils)

    Journal of Computational Physics

    (1974)
  • J. Castillo et al.

    Solution adaptive direct variational grids for fluid flow calculations

    Journal of Computational and Applied Mathematics

    (1996)
  • T.F. Chen et al.

    Numerical construction of optimal adaptive grids in two spatial dimensions

    Computers and Mathematics with Applications

    (2000)
  • J.-C. Yang et al.

    Structured adaptive grid generation

    Applied Mathematics and Computation

    (1994)
  • B. Dacorogna et al.

    On a partial-differential equation involving the Jacobian determinant

    Annales de l’Institut Henri Poincaré. Analyse non Linéaire

    (1990)
  • S. Li et al.

    Moving mesh methods with upwinding schemes for time dependent PDEs

    Journal of Computational Physics

    (1997)
  • L. Chacón et al.

    A 2D high-β Hall MHD implicit nonlinear solver

    Journal of Computational Physics

    (2003)
  • W.J. Rider et al.

    A multigrid Newton–Krylov method for multimaterial equilibrium radiation diffusion

    Journal of Computational Physics

    (1999)
  • D.A. Knoll et al.

    On Newton–Krylov multigrid methods for the incompressible Navier–Stokes equations

    Journal of Computational Physics

    (2000)
  • L. Chacón et al.

    An implicit energy-conservative 2d Fokker–Planck algorithm. II. Jacobian-free Newton–Krylov solver

    Journal of Computational Physics

    (2000)
  • L. Chacón et al.

    An implicit, nonlinear reduced resistive MHD solver

    Journal of Computational Physics

    (2002)
  • D.A. Knoll et al.

    A multilevel iterative field solver for implicit, kinetic, plasma simulation

    Journal of Computational Physics

    (1999)
  • J.U. Brackbill et al.

    FLIP: a method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions

    Journal of Computational Physics

    (1986)
  • J.U. Brackbill

    FLIP MHD: a particle-in-cell method for magnetohydrodynamics

    Journal of Computational Physics

    (1991)
  • D. Sulsky et al.

    A numerical method for suspension flow

    Journal of Computational Physics

    (1991)
  • M. Cho et al.

    r-Adaptive mesh generation for shell finite element analysis

    Journal of Computational Physics

    (2004)
  • W. Huang et al.

    Moving mesh partial differential equations (MMPDES) based on the equidistribution principle

    SIAM J. Numer. Anal.

    (1994)
  • M.J. Baines

    Least squares and approximate equidistribution in multidimensions

    Numerical Methods for Partial Differential Equations

    (1999)
  • J.F. Thompson et al.

    Numerical Grid Generation: Foundations and Applications

    (1985)
  • V.D. Liseikin

    Grid Generation Methods

    (1999)
  • Cited by (0)

    View full text