An optimal robust equidistribution method for two-dimensional grid adaptation based on Monge–Kantorovich optimization
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 be a bounded domain with boundary . We define a two-dimensional coordinate transformation in physical space between the coordinates of an initial grid
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 and the identity mapping between the logical space and the initial grid in physical space: . This implies .
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 .
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 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)
Survey of dynamically-adaptive grids in the numerical solution of partial differential equations
Applied Numerical Mathematics
(1985)Adaptive grid generation
Computer Methods in Applied Mechanics and Engineering
(1987)Equidistribution schemes, Poisson generators, and adaptive grids
Applied Mathematics and Computation
(1987)Variational grid adaptation based on the minimization of local truncation error: Time independent problems
Journal of Computational Physics
(2004)- et al.
Adaptive zoning for singular problems in 2 dimensions
Journal of Computational Physics
(1982) Adaptive grid generation from Harmonic maps
Journal of Computational Physics
(1991)An adaptive grid with directional control
Journal of Computational Physics
(1993)- et al.
An r-adaptive finite element method based upon moving mesh pdes
Journal of Computational Physics
(1999) Variational mesh adaptation: isotropy and equidistribution
Journal of Computational Physics
(2001)- et al.
A fully implicit, nonlinear adaptive grid strategy
Journal of Computational Physics
(2006)