{'fraction_non_alphanumeric': 0.03582350371341197, 'fraction_numerical': 0.011031017911751857, 'mean_word_length': 4.51844262295082, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast Spherical Centroidal Voronoi Mesh Generation: A Lloyd-preconditioned LBFGS Method in Parallel

Centroidal Voronoi tessellation (CVT)-based mesh generation is a very effective technique for creating high-quality Voronoi meshes and their dual Delaunay triangulations that often play a crucial role in applications, including ocean and atmospheric simulations using finite volume schemes. In the next generation climate models, the spacing scales change dramatically across the whole sphere and require ultra-high resolution and smooth transitions from coarse to fine grid regions. Thus fast and robust spherical CVT (SCVT) meshing algorithms become highly desirable. In this paper, we first propose a Lloyd-preconditioned limited-memory BFGS method for constructing SCVTs that is also applicable to the construction of CVTs of general domains. This method is then parallelized based on overlapping domain decomposition, enabling excellent scalability on distributed systems. Results of several computational experiments show that the new method could incur computational time costs one order of magnitude smaller compared with some existing methods for generating large-scale highly variable-resolution meshes, while also providing significantly improvements in mesh quality.


Introduction
A centroidal Voronoi tessellation (CVT) [5] of a given domain is a Voronoi tessellation whose generators coincide with the mass centroids of their corresponding Voronoi cells. Because of its ability to construct point sets that are locally quasi-uniformly distributed with respect to a given point-density function, CVT is a desirable technique for creating high-quality Voronoi and Delaunay meshes [6,7]. For example, such meshes are often crucial for global or regional climate modeling using finite volume schemes [19,20,22,24].
CVTs can be constructed using probabilistic methods such as MacQueen's method [17] or probabilistic Lloyd methods [14] that are elegant random sequential sampling methods that do not require the explicit construction of Voronoi tessellations. However, these probabilistic methods feature very slow convergence so that more efficient deterministic methods for CVT construction are more often used. Typical deterministic CVT construction algorithms include Lloyd's method [16,5], the Lloyd-Newton method [3], and quasi-Newton methods such as limited-memory BFGS (LBFGS) method [15,9].
In addition to the geometric definition given above, CVTs also have an analytical definition as being a critical point of an "energy" functional associated with a set of generators [5]. In particular, a local minimizer of a CVT energy functional determines a stable CVT. Lloyd's method can be viewed as a linearly convergent gradient descent optimization method [5] with respect to the CVT energy function. The Lloyd-Newton method is an attempt to accelerate CVT construction by taking advantage of its quadratic rate of convergence. However, it is not often used for high-resolution CVT grid generation due to the burdensome computation of a Hessian inverse at each iteration.
Quasi-Newton methods are usually more feasible choices because they make use of efficient approximate Hessian matrices or their inverses, and generally have super-linear convergence. In [15] the LBFGS quasi-Newton method and a preconditioned version thereof are used to construct CVTs in planar regions or surfaces. The preconditioner used in [15] is the incomplete Cholesky factorization of a modified Hessian matrix; however, this approach is not guaranteed to be stable or effective [1,9] and, in addition, it has been shown that such LBFGS methods preconditioned in this manner do not apparently improve the efficiency of LBFGS [15] in terms of computational time. A graph Laplacian preconditioner is used in [9] which reduces the computational time to 56.9% on average for constructing CVT meshes with certain non-uniform point-density functions. However, the graph Laplacian preconditioner can be troublesome for generating high-resolution meshes not only because of the limited efficiency due to having to repeatedly solve a large-scale linear system, but also because it is limited effectiveness in reducing the CVT energy.
In recent years, the need for improving the efficiency of the next generation of the computational climate models has driven the need for the more efficient construction of large-scale meshes on the sphere, meshes that feature high-quality and highly variable-resolution with smooth transitions. In this paper, we propose a fast and effective quasi-Newton method for the construction of spherical CVTs (SCVT) but which is also applicable to the construction of CVTs in more general domains. One of our key contributions is the development of a Lloyd-preconditioned LBFGS algorithm; specifically, the Lloyd step is iteratively executed and used as approximations of varying initial Hessian inverses. For quasi-uniform meshes, our Lloyd-preconditioned LBFGS scheme performs similarly to LBFGS. However, for highly variable-resolution meshes, it can dramatically speed up LBFGS while also providing significantly improved mesh quality. Also differing from all previous efforts using LBFGS for computing CVTs, we develop a parallel implementation of our method on distributed systems by using an overlapping domain decomposition approach [12]. Our parallel algorithm features well-balanced loading of mesh points and has excellent performance with respect to strong scaling efficiency. Parallel test results show almost linear speedup when the number of mesh points is of the order of millions.
The remainder of the paper is organized as follows. In Section 2, we review the definitions and basic properties of CVTs and their extension to surfaces/manifolds and also review some popular CVT construction algorithms. Then, in Section 3, we propose the Lloyd-preconditioned LBFGS method for computing SCVTs and its parallelization on distributed systems. Several numerical experiments and comparisons are presented in Section 4 to demonstrate the performance of our method in both serial and parallel environments. Some concluding remarks are made in Section 5.

Review of centroidal Voronoi tessellation and its generalizations
In this section, we first introduce some concepts and properties of centroidal Voronoi tessellation and then review some existing algorithms for their construction.

Definitions and properties of CVTs
Given a bounded domain Ω ⊂ R N and a set of points {z i } K i=1 ⊂ Ω, the Voronoi region V i corresponding to the generator (or site) z i is defined as where · denotes the Euclidean norm in forms a special partition of Ω in the sense that V i ∩ V j = ∅ for i = j and Ω = ∪ i V i . Voronoi cells are convex polytopes except possibly for those whose boundary intersect with the boundary of Ω. The dual graph of a Voronoi tessellation corresponds to the Delaunay triangulation of the set of Voronoi generators. Given a point-density function ρ(x) > 0 defined over Ω, the mass center z * i of a Voronoi cell V i is defined by In general, z i = z * i . In the special case that z i = z * i for i = 1, . . . , K, i.e., each generator coincides with the mass centroid of the corresponding Voronoi cell, the Voronoi tessellation is referred to as a centroidal Voronoi tessellation (CVT) [5].
Alternately and equivalently to the just given geometric characterization of CVTs, one can define them through a variational approach. Letting Z = {z i } K i=1 denote an arbitrary ordered sequence of points in Ω and V = {V i } K i=1 an arbitrary tessellation of Ω, we define the "energy" functional A CVT is a critical point of this functional [5] in the sense that among all possible sets of K points Z in Ω and among all possible tessellations V of Ω into K nonoverlapping subregions, this functional is rendered stationary if and only if {Z, V} form a CVT. Furthermore, local minimizers {Z, V} of F determine stable CVTs. Stable CVTs are more desirable in practical applications than unstable ones which correspond to saddle points of F. Next, consider the tessellation of a compact surface or manifold S ⊂ R N . Given a set of points {z i } K i=1 ⊂ S, we can define the corresponding Voronoi tessellation of S similarly to (1) as For each Voronoi cell V c i , we use a generalized definition for the center of mass to ensure that it lies on the surface S. Specifically, we refer to z c i as the constrained mass centroid [7] of V c i on S if where the point-density function ρ(x) > 0 is defined on the surface S and dσ(·) denotes the surface area element. Note that in general the classical mass centroid z * i defined by (2) does not lie on S. However, as pointed out in [7], there is a connection between the constrained mass centroid z c i and the classical mass centroid z * i , namely, if z c i ∈ S − ∂S, then z * i − z c i is normal to the surface S at z c i , i.e., z c i is the projection of z * i onto S along the normal direction to S at z c i . This feature is especially useful if S is a sphere because on the one hand, the construction of a constrained mass centroid using its definition (4) is considerably more difficult compared to the construction of the classical mass centroid and, on the other hand, for S a sphere, the former is merely the intersection of the radial line going through the latter with the sphere.
The tessellation of S defined by (3) is referred as a constrained centroidal Voronoi tessellation (CCVT) if z i = z c i for i = 1, . . . , K, i.e., the generator associated with each Voronoi cell V c i coincides with the constrained mass centroid of that cell. It is also proved in [7] that a stable CCVT can also be defined as a local minimizer of the CCVT energy function It is worth emphasizing that although the generators of CCVTs are constrained to lie on a surface/manifold, the distance metric used is still the Euclidean distance rather than the geodesic distance.
When S is a sphere in R 3 or a subset of a sphere, as is the case in ocean modeling, a CCVT is referred to as spherical CVT (SCVT).

Some existing deterministic CVT construction algorithms 2.2.1. Lloyd's method
An elegant deterministic method for computing CVTs/CCVTs is Lloyd's method which is a natural byproduct of the geometric characterization of CVTs. Starting from an initial placement of K points in Ω, a Lloyd iteration determines the Voronoi diagram of Ω associated with those points, then moves each point to the mass centroid of its Voronoi cell. This pair of steps is repeated until a sufficient small decrease in the point movement is met. Lloyd's method can be viewed as a fixed-point iteration for the Lloyd map Convergence analyses for some typical cases are given in [5,4].
Based on the variational characterization of a CVT as a critical point of the CVT energy function F(Z, V), a CVT can also be constructed by derivative-based optimization methods because the energy function F is C 2 smooth for a convex domain [15]. Recognizing that in Lloyd's method V = V(Z), i.e., the tessellation V is the Voronoi tessellation corresponding to Z, we have that F = F(Z, V(Z)), i.e., we can view F as a function of Z only. Then the gradient of F(Z, V(Z)) is given by where we use the fact ∂F ∂V = 0 at (Z, V(Z)) due to the property of Voronoi regions [10,5]. Denoting the mass Thus, the Lloyd iteration Z (n+1) = T(Z (n) ) can be rewritten as with the mass matrix M (n) being a diagonal matrix diag(m i ) associated with the Voronoi tessellation V(Z (n) ). In the context of optimization, Lloyd's method is a special case of a gradient descent method. Although it has only a linear rate of convergence, the Lloyd iteration always decreases the energy function F until a local minimizer is obtained without the need of line searches for step-size control.

Newton and quasi-Newton methods
An attempt to accelerate CVT construction is the Lloyd-Newton method [3]. Although the Newton iteration has quadratic rate of convergence, it is usually not feasible for large scale CVT grid generation due to the burden of storing and computation of Hessian inverses. For this reason, quasi-Newton methods are better choices because they use approximations of the Hessian matrices or their inverses and have super-linear rates of convergence. Specifically, given an approximation B (n) of the Hessian matrix H(Z (n) ) of F(Z (n) ), a quasi-Newton iteration operates as follows: for n = 0, 1, 2, . . .
Here q is a search direction along which an optimal step size α (n) is determined by a line search algorithm to satisfy the Wolfe conditions [18]. The choice B (n) = H(Z (n) ) results in the classical Newton method.
One of the most popular quasi-Newton methods is the BFGS method [18]. In the standard BFGS approach, at the n-th step, the inverse of Hessian matrix is iteratively updated by a symmetric matrix H (n+1) being close to the last Hessian inverse approximation and satisfying the secant equation Such derivation forces the quadratic forms of function F around Z (n) and Z (n+1) match gradients respectively at these two steps. The BFGS updating formula in the product form can be given explicitly as H (n) = S T n−1 · · · S T n−m(n) H (0) S n−m(n) · · · S n−1 + S T n−1 · · · S T n−m(n)+1 ρ n−m(n) s n−m(n) s T n−m(n) S n−m(n)+1 · · · S n−1 + S T n−1 · · · S T n−m(n)+2 ρ n−m(n)+1 s n−m(n)+1 s T n−m(n)+1 S n−m(n)+2 · · · S n−1 + · · · + ρ n−1 s n−1 s T n−1 , In the classic BFGS updating scheme, m(n) = n is taken and the initial approximation H (0) is often set to be a constant multiple of the identity matrix.
In the limited-memory BFGS (LBFGS), one needs to specify M , the number of BFGS corrections needed to be stored. That is, only the recent M pairs of {s i , y i } are kept. The BFGS approximation is modified by setting m(n) = min{n, M } and allowing the initial approximation H (0) to vary from iteration to iteration. For instance, H (0) is typically replaced by H is close to the reciprocal of an eigenvalue of H(Z (n) ). In practical computation, the product H (n) ∇F(Z (n) ), rather than the matrix H (n) , is recursively updated. Specifically, the LBFGS updating at step n is done as follows:

Lloyd-preconditioned LBFGS in parallel
In this section, we will first develop a Lloyd preconditioner for the LBFGS method for SCVT construction so that its performance can be further improved, after which we propose an efficient parallelization algorithm for the method on distributed systems based on overlapping domain decomposition, which is important in terms of both computational and storage efficiency.

Lloyd-preconditioned LBFGS
The idea of a preconditioned LBFGS method (P-LBFGS) proposed in [13] is to choose a matrix B (n) very alike to the Hessian H(Z (n) ) and periodically replace line 4 in the LBFGS iteration by It requires that B (n) should be inexpensively achievable, and most importantly, the above equation can be easily solved. If the matrix B (n) can capture significant features of the Hessian, then the convergence of LBFGS can be hastened. It has been shown that the CVT energy function has second-order smoothness for convex domains [15] and that the Hessian matrix can be computed exactly; see [10,15]. Denote by J i the set of node indices of the Voronoi cells adjacent to V i (the Voronoi cell corresponding to generator z i ). Let the coordinates of the point z i (resp. y) be denoted by z ik (resp. y k ), k = 1, . . . , N . Then, the second-order derivatives of the CVT energy are given by the following explicit formulas: for k, l = 1, . . . , N , In [15], a modified Cholesky factorization of the Hessian matrix H(Z (n) ) is applied to obtain the positive-definite matrix B (n) , and perform (6) every N iterative steps, where N is a user-defined parameter. The numerical experiments in [15] show that LBFGS and P-LBFGS yield better performance compared to Lloyd's method and Newton's method. However, the preconditioned version of LBFGS used in [15] only improves LBFGS a little, and these two perform very similarly even in the case of non-uniform point-density functions.
Another approach at developing a P-LBFGS method is to use a graph Laplacian preconditioner (graph Laplacian P-LBFGS), proposed in [9]. They consider CVT construction in a polygonal domain Ω in two dimensions. Let T ij denote the triangle inside V i containing the point z i and the edge V i ∩ V j . If V i has an edge e ib in ∂Ω, denote the triangle formed by e ib and z i as T ib . Then, the graph Laplacian matrix A = (a ij ) is given by According to the numerical tests in [9], for a constant point-density function, the graph Laplacian P-LBFGS has similar efficiency compared to the classical LBFGS algorithm, but for non-uniform point-density functions, the graph Laplacian preconditioner can reduce the computational time by 56.9% on average. However, the use of graph Laplacian preconditioners to speed up LBFGS could be limited for ocean grids due to some difficulties: first, solving the Laplacian linear system on large-scale meshes at each iteration is computationally expensive; second, the Laplacian matrix A is singular for closed surfaces such as the surface of the sphere and thus special care is needed in choosing an appropriate solution for the preconditioning process. Furthermore, the graph Laplacian preconditioner sometimes does not work well in significantly reducing the CVT energy, as we reported in Section 4.2.
In this paper, we propose a new preconditioner (initial Hessian inverse, line 4 in the LBFGS iteration), inspired by Lloyd's iteration, as We call this the Lloyd-preconditioned LBFGS (Lloyd P-LBFGS) method because this approach is equivalent to taking the initial BFGS direction to be the Lloyd's updating direction (see (5)) and the initial line-search step size as 1.
There are several reasons behind this idea. First of all, because the Lloyd map T is continuous, the Lloyd iteration converges globally to a critical point of F if the iterations stay in a compact set [4]. Although the compactness of the iteration has not been rigorously justified in the literature, it seems to be intuitively true and matches practical experiences. This nice property of Lloyd's method is not shared by Newton or quasi-Newton methods. Second, the initial Hessian inverse given by (8) is very simple while still capturing a lot of Hessian matrix information. Note that 2M (n) is the main diagonal part of the Hessian matrix as shown in (7). A connection of (8) to the graph Laplacian matrix is more obvious: it is an approximation of the diagonal components of the Laplacian matrix A because and we can treat . Lastly, the diagonal matrix M (n) itself contains significant information about the density distribution.
To apply LBFGS for minimizing the energy function F c in the case of CCVT, in [15] a "pseudo" derivative is taken by projecting the gradient of the CVT energy onto the tangent plane of the compact surface S. That is, they proposed to use the tangential derivative in the minimization process, where N(z i ) is the normal to the surface S at z i . The updated set of generators needs to be projected back onto the surface at each iteration. Following this approach, we can derive the Lloyd preconditioner for the case of the unit sphere analogous to that for planar domains. Note that Setting c i = V c i yρ(y)dσ, the tangential derivative can be written as The updating by d i is then a Lloyd iteration because the projection of d i onto the unit surface coincides with that of c i . As such, we can take the Lloyd preconditioner in R K×K as

Parallel implementation
Our parallel implementation of the Lloyd-preconditioned LBFGS method for SCVT mesh generation is based on overlapping domain decomposition.

Parallelization of Delaunay triangulations on sphere
As is true for all other deterministic methods for CVT/SCVT construction, the Lloyd-preconditioned LBFGS algorithm involves explicit integrations over Voronoi cells. The Voronoi cells can be determined from the dual Delaunay triangulation. In this subsection we present a parallel algorithm for computing Delaunay triangulations and spherical Delaunay triangulations first studied in [11,12].

Planar Delaunay triangulations. A Delaunay triangulation for a set of points {z
is strictly inside the circumcircle of any triangle in T . Delaunay triangulations maximize the minimum angle of all the angles of the triangles in the triangulation; they tend to avoid extremely acute angles of the triangles. The uniqueness of Delaunay triangulation for {z i } K i=1 is guaranteed if the set of points is in general position (no set of four points lie on a circle whose interior does not contain points from ). The parallelization of planar Delaunay triangulation algorithms has been widely studied in the past decades; see, e.g., [2,8]. Popular algorithms commonly split the global point set into equally distributed subsets, each of which can then be triangulated simultaneously in parallel. The resulting local triangulations are then stitched together to form a global triangulation. Such a merging step may need to modify significant portions of the local triangulations, thus it is the most significant step that affects time efficiency. In this work, we resort to an alternative merging step that needs no modifications of the local triangulations and thus performs well in parallel. The algorithm for planar Delaunay triangulation consists of the following three steps (see Figure 1-left).

Given a planar domain Ω and a set of points {z
in Ω, we preform an overlapping covering of Ω and record the overlapping connectivity or neighboring list. A simple example of the covering C can be a set of p disks denoted by {D l (p l , r l )} p l=1 , where p l and r l are the disk center and radius respectively. The K points in Ω are then distributed into p subregions, each of which is handled by a single processor. Because of the overlapping feature of the covering, some points may belong to multiple processors.
2. On the l-th processor, a local Delaunay triangulation T l is carried out independently. To this end, one can employ a serial Delaunay algorithm such as the Delaunay triangulator in the Triangle software package [23].
3. On each processor, we only keep triangles whose circumcircles are completely contained inside the corresponding subregion. Denoting the new triangulation on the l-th processor by T l , then the circumcircle of any triangle in ∪ l T l does not contain any points from {z i } K i=1 . Each local triangulation T l is then exactly a portion of the global Delaunay triangulation due to the uniqueness of Delaunay triangulations.
Spherical Delaunay triangulation. We now consider spherical Delaunay triangulation on the unit sphere S in R 3 . In this case, triangles and Voronoi cells are replaced by spherical ones whose boundaries are geodesic arcs. The planar parallelization scheme is combined with stereographic projections. A stereographic projection is a conformal mapping that projects a sphere onto a plane tangent to the sphere. Letting t denote the contact point of the sphere and its tangent plane, the stereographic projection of a point z on the unit sphere to the point x on the tangent plane is given by (in Cartesian coordinates) .
The conformality of stereographic projections implies the preservation of circumcircles of triangles along with their interiors and hence the Delaunay criteria. Therefore, the Delaunay triangulation of a subregion of the sphere can be obtained by the stereographic mapping from the planar Delaunay tessellation on the corresponding tangent plane. The adaptation to spherical Delaunay triangulation from a planar one is described by the following three steps; see Figure 1.
1. We preform an overlapping covering of the unit sphere S and record the overlapping connectivity or neighboring list. The covering U could simply be a set of geodesic disks {U l (t l , r l )} p l=1 , where t l is the disk center and r l is the geodesic radius. Points {z i } K i=1 on S are then distributed into p subregions with overlaps.
2. On the l-th processor, we project the point set {z i } K i=1 ∩ U l (t l , r l ) onto the plane tangent to the sphere at t l . We then proceed to construct a planar Delaunay triangulation T p l . The local spherical Delaunay triangulation T s l is then produced by stereographic projection from T p l .
3. In each subregion, say U l (t l , r l ), we only keep spherical triangles whose circumcircles are completely contained inside, so as to guarantee the Delaunay property with respect to the global point set {z i } K i=1 . In particular, for disk subregions, let p c j and r c j denote the center and radius of the circumcircle of the j-th spherical triangle in U l (t l , r l ); the triangle selection criteria is arccos(t l · p c j ) + r c j ≤ r l . Each local spherical triangulation is now exactly a portion of the global spherical Delaunay triangulation.
To obtain optimal load balancing and robustness, the target domain should be decomposed based on the given point-density function. A precomputed ultra-coarse CVT is then a natural choice for effecting the partition. The main challenge is now the determination of the sizes of the overlaps. The overlapping regions should be large enough to make sure no true Delaunay triangles are missing in the final triangulation. A simple overlapping decomposition for quasi-uniform meshes is to use disk regions centered at partition cell centers with radius being the maximum distance from the center to its adjacent centers. More challenges exist for variable resolution meshes. In this work, we take the Voronoi sort method, i.e., the union of the owned partition cell and its immediately adjacent partition cells defines the overlapping domain decomposition for parallel triangulation. Failure may occur when the global point set does not include enough points for the decomposition of choice. A practical suggestion on the number of grid points is at least 16 multiples of the number of partitions, or a number corresponding to at least two-level bisections from the partition CVT, as pointed out in [12].

Parallelization of Lloyd P-LBFGS
To parallelize the Lloyd P-LBFGS scheme, we first transform the global point vector Z into non-overlapping local vectors, so that arithmetic operations like additions and dot products can be efficiently applied on the distributed vectors. This non-overlapping decomposition can be initially carried out according to the disjoint partition cells C 1 , . . . , C p ; see Figure 2-left. On each processor, the global IDs of points in the distributed vector are recorded and will be used for vector assembly after point updating. For any point, say z i , if it is initially in partition cell l by disjoint partition, we call cell l its arithmetic position. As the iterations proceed, z i could move from one partition cell to another, and we call the current partition cell it belongs to as its geometric position. Whereas the geometric position could vary, we always fix the arithmetic position for convenience of computations inside LBFGS. The evaluations of the energy function, the gradient, and the Lloyd preconditioner require integrations over Voronoi cells. All integrations are actually computed based on the dual Delaunay triangulation. Given a triangle (z a , z b , z c ) in the Delaunay tessellation, with its circumcenter c and three edge midpoints z ab , z bc , and z ca , the triangle can be split into six sub-triangles, each made of an original triangle vertex, an edge midpoint, and the circumcenter. The sub-triangle (z a , z ab , c), for instance, is part of the Voronoi cell whose generator is located at z a . Should the orientation of sub-triangles be taken into account, an integration over a Voronoi cell is the sum of integration over such sub-triangles inside the cell. The evaluation of the functional gradient and the Lloyd preconditioner described above are performed in parallel and stored according to the geometric positions of grid points. To be clear, we will call the partition cell associated with which current processor is handling as "my cell", and the union with its overlapped layers as "my region". After points are assigned into "my region", local Delaunay triangulation is performed. We then integrate over "my cell" for energy function evaluation. In the meanwhile, at each point in "my cell", we compute the function gradient and the Lloyd preconditioner. Non-blocking communication is then used to assemble their arithmetic vectors, whose distribution is the same as the non-overlapping partition of the initial global points.
It is worth noting that in our code there is no need for global gathering of local Delaunay triangulations or updated mesh points (Voronoi generators). The communications for updating mesh points are only among two-level neighboring partition cells, under the assumption that grid points move at most into their neighboring partition cells. Here the "two-level" approach is to ensure the updating of points in overlapped layers. The assumption is generally fulfilled in practice by consecutively performing optimization-bisection until the final level of optimization; this is detailed in next section. As such, the communications for assembling the gradient and preconditioner are only in one level, i.e., only among neighboring partition cells. When the assumption is not satisfied, a reset on P-LBFGS is executed (see the illustration in Figure 2-right), i.e., the arithmetic position for each mesh point is re-computed.

Description of the parallel algorithm
There are still several ingredients that need to be considered for parallel P-LBFGS for SCVT construction, such as what are the point-density function, the initialization, and the stopping criteria.
The point-density function ρ(x) plays a crucial role in how the converged mesh points are distributed. For the numerical tests in this paper, we use analytic point-density functions that feature significant physically relevant information drawn from climate modeling. A point-density function of course also can be formulated from a given initial mesh, or even from a local mesh of a subregion of the sphere.
The convergence of iterative methods is often highly sensitive to the initial configuration. Exploring a good initialization for SCVT computation has been an interesting topic in the field. To the best of our knowledge, a Monte Carlo initialization based on the point-density function is often suitable for coarse mesh generation. For high-resolution meshes, however, a consecutive optimizationbisection procedure is almost mandatory to obtain an initial guess that is good enough to save computational cost.
The stopping criteria for P-LBFGS here are based on the maximum number of iterations, the norm of gradient of the energy function, and the step size of mesh point movement.
In summary, the parallel version of the Lloyd P-LBFGS algorithm for SCVT construction consists of the following steps.
1: Build the partition information from a coarse CVT; 2: Set n ← 0; 3: Decompose Z (n) into disjoint partition cells accordingly; store global IDs of points; 4: while stopping criteria not satisfied do 5: Preprocess: assign points from Z (n) into "my region" (2-level neighboring comm. ); go to Step 3 if a point in Z (n) moves out of its initial neighboring partition cells for points in "my region", build local Delaunay triangulation; 6: For points in "my cell", compute F(Z (n) ), ∇F(Z (n) ), H (0) n ; 7: Sum F(Z (n) ); assemble ∇F(Z (n) ), H (0) n to arithmetic positions (1-level neighboring comm. ) 8: Compute search direction q by P-LBFGS; 9: Update Z (n+1) = Z (n) + α (n) q by line search: occasionally perform Step 5 on Z (n) + sq and evaluate F, for some step size s.

Performance tests
Our code for the Lloyd-preconditioned LBFGS solver for SCVT construction is written in C++ and MPI; the Delaunay triangulation is handled by the "Triangle" package, the LBFGS routine is from the serial HLBFGS package by Yang Liu (in which the line-search routine is converted from the original LBFGS Fortran code by Jorge Nocedal). We set the maximum number of function evaluations for the line search to ten and store at most seven BFGS corrections in memory (M = 7), as suggested in [15]. The SCVT grid generator in this study is under active development and the code is available at https://github.com/hyang52/scvt HLBFGS.
The quality of the Voronoi cell V of the generated meshes is quantified by cellQ cellQ(V ) = min cell edge max cell edge and the triangle quality triQ, which is defined on an triangle with edge lengths a, b, c by Three different point-density functions defined on the unit sphere are used in our tests to represent some typical cases. For each point-density function, we use the max-min ratio of the pointdensity to control the max-min ratio of mesh size because it is conjectured that for any two Voronoi cells V i and V j from a SCVT, their mesh spacings are related to the point-density function as The first point-density function we used is given by , so that the mesh size on the equator is three times of that on the poles. This size distribution is typically used in some simulations of ocean modeling. We will call it the X3 point-density; see the SCVT mesh in the left column of Figure 3. The other two point-density functions are defined as where d(z i ) is a distance function that can be simply the geodesic distance from z i to a pre-specified fixed point z c on the sphere. This definition results in a relatively large value of ρ within a distance β from the point z c . The point-density transitions from a large value to small value across a geodesic distance of α.
To have a variety of tested densities, we define d(z i ) in different ways for the second and third point-density functions. Specifically, for the second point-density function, we define where z ic is a point whose latitude and longitude coincide with z i and z c , respectively, and vice versa for z ci . The norm | · | S denotes geodesic distance. We choose the widths w lon and w lat to be 0.3 and 1.2 with z c = (1.0, 0.0, 0.0), and set β = π/6 and α = 0.3. The minimum point-density is set to be γ = (1/16) 4 , so we refer to the defined function as the X16 point-density; see the SCVT grid in the middle column of Figure 3). For the third point-density, we take a function used in [19] for a multi-resolution modeling with the shallow water equations. The distance is now d(z i ) = |z i − z c | S with z c = (0.0, −0.866, 0.5), the parameters are taken as β = π/6 and α = 0.15. Here, we set the minimum point-density γ = (1/64) 4 and refer to the function as the X64 point-density; see the SCVT grid in the right column of Figure 3). There are two reasons of using this X64 point-density. On the one hand, it corresponds to a highly variable multi-resolution mesh that is needed in global ocean and coastal system modeling for which the mesh size may range from a O(10 4 )km global scale to O(10 −1 )km local scale. On the other hand, it can be used in global spherical mesh generation for the purpose for limited-area climate modeling, if one culls the coarse portion from the global one.

Comparison of Lloyd, LBFGS and Lloyd P-LBFGS in serial
Lloyd's method is the best-known algorithm for SCVT construction and LBFGS is a sound way to speed up the Lloyd iteration, but we believe that the Lloyd P-LBFGS method will further speed up LBFGS, especially for highly variable multi-resolution meshes. To compare the performance of these three methods, we execute a sequence of tests for different scenarios for serial implementations. In these tests, we only work on coarse meshes with 2562 and 10242 SCVT generators and all the runs are effected on a PC with a 2.70GHz Intel i7-3740QM CPU. The mesh points are initialized by Monte Carlo sampling based on the corresponding point-density functions. To compare both the efficiency and the mesh quality, we set the maximum number of iterations to be 2,000 and the tolerance on the size of the point movement to be 5 × 10 −4 . Except for that, the quasi-Newton iterations will stop when ∇F n /F n ≤ 5 × 10 −4 or |F n − F n−1 |/F n−1 < 10 −7 . The tolerance 10 −7 on the control of the energy decrease is to avoid wandering. For the integration over a triangle, we use a 4-point quadrature rule in the X3 and X16 cases, and a 9-point quadrature in the X64 case which includes extremely coarse cells. Table 1 provides detailed information of the three tested methods with 2,562 generators. Because the SCVT energy functional has many local minimizers, different optimization methods with even the same initial guess can still lead to different minimizers. For this reason, the final gradient norm ∇F plays a role in judging the quality of final SCVT mesh.  Table 1 shows that for all three point-density functions, Lloyd P-LBFGS is much more efficient than LBFGS and Lloyd (1.18 to 5.02 times faster than LBFGS and 5.66 to 7.19 times faster than Lloyd). We also observed that the relative efficiency of Lloyd P-LBFGS becomes even better in the case with 10,242 generators, as illustrated in Table 2. Among the three methods Lloyd P-LBFGS always returns the smallest value on the energy function and nearly so for the gradient norm, which usually mean better mesh quality. The speed-up of LBFGS by the Lloyd preconditioner is more obvious as we increase the mesh size ratio; this is also shown by the plots of the iteration histories in Figure 4 which provides the numerical results in terms of the energy functional and the gradient norm with respect to the number of iterations and the computation time.
SCVT meshes with 2,562 generators generated by Lloyd P-LBFGS are plotted in Figure 3 from which we can see that the mesh points are very regularly distributed and that the mesh transitions from coarse to fine regions are very smooth. The quality of SCVT meshes can also be characterized by the triangle quality measure ranging from 0 to 1; see the definition in (9). In Figure 5, we plot the distribution of the triangle quality of the SCVT meshes with 10,242 generators produced by the three methods; we here only plot the case of X3 and X16 with 10,242 generators because the X64 case has some very triangular cells that are too coarse to be properly visualized by the Paraview software. In this figure, the lighter the color, the better the quality. The results show that in general both LBFGS and Lloyd P-LBFGS produce better SCVT meshes in term of quality than Lloyd and for the X16 case, the Lloyd P-LBFGS clearly provides meshes with better quality than other two

Comparison of Lloyd and graph Laplacian preconditioners in P-LBFGS in serial
We next make a direct performance comparison of the Lloyd preconditioner with the graph Laplacian preconditioner in the P-LBFGS scheme. The graph Laplacian is computed with a Cholesky decomposition solver as done in [9]. Because the Laplacian matrix is singular, the numerical solution is unique up to a constant; depending on which graph Laplacian preconditioner is used, the corresponding P-LBFGS method could perform very differently. Some solutions of the Laplacian system may lead to a direction with a poor length for line search or even a non-decent direction. As a result, we will handle the singularity of the Laplacian operators by a diagonal perturbation. Two approaches are tested . One is the diagonal shifting by adding scaled identity: A ← A + 10 −6 I. We call this approach Laplacian(a6). The other one is to multiply the diagonal entries by a constant factor: A ii ← A ii · (1 + 10 −2 ). This method is denoted by Laplacian(m2). We apply the graph Laplacian preconditioner in each iteration so that the numerical results can represent, in principle, the best performance of Laplacian preconditioner.
It is worth mentioning that the point-density function in [9,Example 4] has min-max ratio e −80 , thus the mesh size ratio is as huge as e 20 ; similarly in [9,Example 5], the mesh size ratio is e 2.5 ≈ 12.18. Our testing densities X16 and X64 are in a category comparable to [9]. In the following tests, the parameter configuration is exactly the same as in Section 4.1. We investigate, in different scenarios, not only the necessary number of iterations for convergence but also the final mesh quality which is an even higher priority in many applications [19,21].
Detailed test results are provided in Table 3 and Figure 6. We see that Laplacian(m2) has the overall fastest convergence rate if we only focus on the number of iterations without mentioning the computational time. Laplacian(a6) requires fewer iteration than the Lloyd preconditioner in the X3 test case, whereas not significantly fewer in the X16 case, and, in the X64 test case, it requires more iterations. Thus, we take Laplacian(m2) as a better representative for the graph Laplacian preconditioners proposed in [9]. The graph Laplacian preconditioner can dramatically reduce the necessary number of iterations in P-LBFGS; this is their main advantage as it was proposed. A natural concern is the computational cost incurred by having to solve discrete Laplacian systems. We observe that with respect to computational time, the Lloyd    in all test cases. Another aspect which we care even more about is mesh quality. As one can see from Table 3 and Figure 6, Laplacian(m2) always returns larger values than Lloyd P-LBFGS for the final energy and the gradient norm. We also find that in this set of tests Laplacian(m2) actually always stops based on the criterion of |F n − F n−1 |/F n−1 rather than on the gradient criterion. It seems the final energies only have little differences, but the mesh qualities are significantly different, as illustrated by the plot of cell quality in the third row of Figure 6. We conclude that Lloyd preconditioner outperforms the graph Laplacian(m2) at providing a SCVT grid with "optimal" quality and computational times.

Parallel scalability of the Lloyd P-LBFGS
Now we study the scaling performance of the parallel P-LBFGS method. All the parallel tests were run on a high performance computing (HPC) cluster with 16 Intel Xeon "model E5-2670" cores per node.
We first present in Figure 7 the speedup of parallelization with varying number of generators (163,842, 655,362, 2,621,442) and different point-density functions. The speedup is measured    scalability ends up being sub-linear. As the number of generators increases, the speedup gets better and better. When the number of grid points is in million scale, the parallel solver features almost linear scaling up to 256 processors. We next show the scaling results towards generating the SCVTs with 655,362 generators. Initialized from a convergent CVT mesh of 2,562 generators, our code consecutively bisects and optimizes the mesh until a converged SCVT with 655,362 generators is obtained. The iteration history is recorded for the 655,362 generator mesh and also the 40,962 generator mesh generated along the consecutive bisection procedure. To measure the parallel efficiency, we launch the program with different number of processors ranging from 1 to 128. Figure 8 shows the log-plot of the running time t against the iteration number n. As they are related by t = a · n + c, when n is relatively large we can approximate log t ≈ log a + log n. Thus, the distance between the right ends of adjacent curves reflects the gain on time efficiency per iteration when the number of processors is doubled. Excellent scalability is confirmed because those distances are fairly evenly distributed. The plots of F and ∇F vs. the running time for the 40,962 points mesh are shown in Figure 9.  Finally, in Figure 10, we plot the distribution of triangle qualities of the SCVT meshes with 655,362 generators created by the parallel Lloyd P-LBFGS method with 128 processors. The total running times are respectively 1296.4 seconds for X3 point-density, 1939.6 seconds for X16 pointdensity, and 1914.9 seconds for X64 point-density.

Summary
A Lloyd-preconditioned LBFGS method and its parallel implementation are proposed and numerically investigated for SCVT computation. Whereas LBFGS often speeds up the CVT computation as compared to the classical Lloyd method, it may lose "optimal efficiency" in generating highly variable multi-resolution meshes. The Lloyd-preconditioned LBFGS method can further speeds up the CVT computation and overcome this difficulty in large scale SCVT computation.
For quasi-uniform meshes, the Lloyd-preconditioned LBFGS method performs similarly as LBFGS, which is also the case when using the graph Laplacian preconditioner. For typical multi-resolution 1-proc 2-proc 4-proc 8-proc 16-proc 32-proc 64-proc 128-proc  meshes, however, the Lloyd-preconditioned LBFGS method significantly improves the performance of LBFGS, while also producing better mesh quality. Because the method enables efficient SCVT generation in highly variable multi-resolution, it is also useful for limited-area climate modeling by culling coarse grids. Note that it is often challenging to maintain regularity of meshes at boundaries by direct local mesh generation. The parallel Lloyd-preconditioned LBFGS algorithm features well-balanced loading of grid points and performs well in the parallel scalability tests. It enables more convenient (faster mesh generation) and more stable (better mesh quality) global climate modeling, especially for meshes with ultra high-resolution. The proposed efficient CVT grid generator will thus play a fundamental role in next generations of climate modeling, where the spacing scale could be as small as 100 meters. Our interests will also include more accurate numerical simulations in ocean and atmosphere modeling based on SCVTs.