A gradient augmented level set method for unstructured grids

Journal of Computational Physics 258 (2014) 47–72

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

A gradient augmented level set method for unstructured grids Arne Bøckmann ∗ , Magnus Vartdal University of Oslo (UiO), Department of Mathematics, PO Box 1053, Blindern, NO-0316 Oslo, Norway

a r t i c l e

i n f o

Article history: Received 4 October 2012 Received in revised form 26 June 2013 Accepted 14 October 2013 Available online 21 October 2013 Keywords: Level set Unstructured grid Compact stencil Gradient augmented Fast marching method

a b s t r a c t We present and test a gradient augmented level set method for unstructured grids, which has been designed to allow for easy integration in parallel flow solvers. The means to this end is a Hermite formulation, where the gradient is stored and advected as an independent variable, along with the level set function. The method uses narrow-band storage in conjunction with a Fast Marching Method, which is also designed to take advantage of the gradient augmentation. We demonstrate through standard passive advection test cases that the proposed method is able to compute interface motion with performance comparable to popular high-performing methods on structured grids. © 2013 Elsevier Inc. All rights reserved.

1. Introduction Modern fluid simulation of practical problems usually involve complex geometries, variable spatial resolution and parallel computing. The introduction of phase interfaces adds to the complexity. Whereas there exists an abundance of methods and algorithms to track fluid interfaces, much fewer possess the generality that allow them to be implemented with an unstructured, parallel solver. Volume Of Fluid (VOF) methods and Level Set methods are two of the most popular interface capturing approaches used on stationary grids. In the original VOF method [1], the volume fraction of fluid in the computational cell is used to identify and track the surface, ideally conserving mass. However, VOF encounters a common problem: a characteristic of a fluid surface or interface is that it represents a discontinuity, which does not readily lend itself to polynomial interpolation schemes. While several variants of shock capturing schemes exist, these introduce significant numerical diffusion when applied to discontinuities. Thus, in [1], the surface is reconstructed in a geometric fashion, where chunks of fluid are transfered from one cell to the other. To maintain a sharp interface, additional ad-hoc procedures are introduced after advection. This is a relatively crude process, and it is particularly challenging on unstructured grids. Thus, while conserving mass, VOF has in many cases difficulties with accurately representing the motion of an interface; see [2] for a modern example. Level set methods, on the other hand, convert the scalar function, in which the interface is embedded, to a form more suitable for polynomial interpolation schemes. This is done through the solution of a reinitialization equation – most often an eikonal equation – which converts the level set function into a signed distance function. Since the distance to the interface is not a conserved quantity, level set methods using this approach do not conserve mass exactly. While the signed distance function represents the interface in a relatively smooth manner, numerical diffusion is still present, and this often causes unacceptably high local mass loss or gain, especially in regions of high curvature. Several strategies exist to counteract this. The conceptually simplest is to increase the resolution of the grid near the interface, either through adaptive mesh refinement [3], or a finer, separate background grid [4]. Another approach is to

*

Corresponding author. E-mail address: [email protected] (A. Bøckmann).

0021-9991/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.10.024

48

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

increase the numerical accuracy. The WENO-type high-order, shock-capturing schemes [5] are popular for structured grids. Implementation on unstructured grids is however more complicated [6], and a general drawback of WENO schemes is that the numerical stencils are rather large, thus complicating data transfer in parallel settings. A simpler approach is the use of Lagrangian particles to track the interface. The Particle Level Set method [7] is one such method, where particles of finite radius are placed tangent to the interface to correct numerical error in the Eulerian solution. By using Lagrangian interface tracking in conjunction with the level set method, topology changes may be simplified as compared to a purely Lagrangian method. While the original Particle Level Set method used a WENO scheme for the grid-based advection, it was discovered in [8] that since the accuracy mainly depends on the particle correction, a first order semi-Lagrangian scheme for the grid-based quantities is sufficient. Lagrangian augmentation of the interface advection has many attractive features. But like all of the discussed methods, it introduces an additional layer of complexity. Moreover, as the interface is stretched or compressed, so are the particles. This necessitates particle seeding or deletion. When seeding new particles, some kind of interpolant must be used to approximate the interface position from the grid data, and the accuracy of the interface reconstruction is limited by this interpolant. Due to the relatively limited information contained in each particle (position, radius and sign), and the fact that the corrections are made on a per-particle basis, a relatively dense particle distribution is required. In [7], a default of 64 particles per cell per phase was used in 3D. Particles were seeded in a band that spanned around three cells from the interface, yielding a large number of particles. In [9], gradient and curvature information were also embedded in Lagrangian particles and used in conjunction with a level set method. The information transfer from particle to grid was however done by a first order Taylor expansion only, and it remains unresolved how reseeding and topology change should be treated in a robust manner. If local mass conservation is a pressing issue, the signed distance function may be replaced by a smeared-out approximation of the volume fraction [10,11], however, this may decrease the accuracy of the numerical schemes. Global mass conservation can always be assured by iterating a scheme that moves the interface in the normal direction [12]. This may be an acceptable approach if the local mass loss/gain is low, but over longer periods of time leads to unphysical results. As mentioned above, high-order schemes can reduce the mass loss, but these often require large stencils. Cubic Hermite interpolation is one of the exceptions – all data required to construct the interpolant is contained within a single cell. The gradients may be obtained by a transport equation, which is derived by taking the gradient of the scalar advection equation, or its discrete analogue. While maintaining compact stencils, this also introduces additional effective resolution, at the expense of increased data storage. Transporting the gradient as an independent quantity has been used in Cubic Interpolated Pseudo-particle/Constrained Interpolation Profile (CIP) [13,14] methods for some time, however, only recently has this approach spread to level set [15] and VOF methods [16]. In [15], a semi-Lagrangian scheme based on tensor product cubic Hermite splines was used. This required a structured grid, and the issue of reinitialization was not treated. Specifically designed for use with the advection scheme proposed in [15], a reinitialization procedure using a Fast Marching Method (FMM) was introduced in [17]. Combined, these two should constitute a good choice for structured grids. In [18], flow simulations using semi-Lagrangian gradient transport and the tetrahedral cubic Hermite interpolant described in [19] were presented, though interface advection was not treated. Other interesting approaches to high-order semi-Lagrangian advection schemes are found in [20,21]. In [20], a method for structured grids which uses a narrow-band formulation with arbitrary order Gauss–Lobatto stencils is presented. In [21], a conservative, 2D unstructured scheme is formulated using a combination of point values updated by a semi-Lagrangian scheme, and area average updated by a Finite Volume scheme. These are only a few examples of the versatility of semi-Lagrangian schemes. An additional advantage is that they are not limited by the Courant–Friedrichs–Lewy (CFL) condition, as the characteristics may be traced back across several cells during one time step [22]. In the present work, we extend the work of [15] to general 3D unstructured grids, and include reinitialization of the signed distance property, which is integrated into the gradient-augmented setting. This is done by a Fast Marching Method which computes both the level set function and its gradient. For the FMM, we apply Huygens’ principle to obtain a local approximation of the distance function. This approach requires local solutions of a small optimization problem, and differs from the Finite Difference (FD) approach used in [17]. Huygens’ principle leads to the same type of update rule as presented in e.g. [23,24], which are based on a control-theoretic approach. We use a narrow-band method with storage based on a tetrahedral subgrid, meaning that general, unstructured primal grids are decomposed into tetrahedra. To avoid overshoots and generation of artifacts in high-curvature areas, we develop a gradient limiting algorithm which is demonstrated to be important for the overall accuracy. The proposed method can readily be incorporated into a parallel setting, as we adhere strictly to compact stencils. While this has been done, the detailed parallel implementation and performance are not the focus of the present work. It should be noted that in the Finite Element framework, another approach with many of the same attractive features we seek already exists. The Discontinuous Galerkin Method has been applied with the level set approach on unstructured grids [25,26] with excellent results. The rest of the article is organized as follows. A brief description of the level set method, narrow-band and grid structure is given in Section 2. Cubic Hermite interpolants for 1-, 2- and 3-simplices are described in Section 3, followed by an explanation of the semi-Lagrangian advection scheme for function value and gradient in Section 4 and the reinitialization procedure for the level set function in Section 5. We test the method on standard passive advection cases in Section 6, and discuss the overall method and results in Section 7. To promote readability, much of the discretization details, examples and theoretical support are left to Appendices.

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

49

Fig. 1. Grid structure.

2. Level set method and grid structure The level set method, introduced in [27], treats interface advection in n dimensions by embedding the (n − 1)dimensional interface in an n-dimensional function Φ(x, t ), here referred to as the level set function. The interface is described by an isosurface of Φ . This dictates that the material derivative of Φ equals zero, or in other terms, that Φ is advected passively:

∂Φ + u · ∇Φ = 0, ∂t

(1)

where u is the advecting velocity field. Using the zero-isosurface, the set of points Γ describing the interface are given by

  Γ (t ) = x ∈ Ω : Φ(x, t ) = 0 ,

(2)

where Ω is a domain contained in R . Since the zero-isosurface of Φ is the main interest, values of Φ away from the zeroisosurface play a less important role. In practical, discretized applications, Φ can be chosen freely away from the interface as long as it has a diminutive or no effect on the zero-isosurface found by interpolation, and consequently as a minimal requirement has the correct sign. As values of Φ outside the extent of the interpolation schemes used to describe the zeroisosurface do not affect the zero-isosurface, it suffices to advance Eq. (1) in a narrow band across the zero-isosurface [28]. This can greatly reduce the computational cost and storage requirements. To prevent numerical diffusion from propagating large irregularities from new members of the band, the linear Fast Marching Method described in Section 5 is applied to assign monotonically increasing values of |Φ| to the new band members, based on values in the existing band members. Instability near the edge of the band was also noted as a potential problem in [29], which used smoothing to prevent it. A cut-plane of prismatic triangular cells containing an interface in a narrow-band structure is shown in Fig. 1(b). The band is updated after every advection step; new members are added, and existing members are removed to maintain a constant thickness in terms of cells. The storage points in the band are the grid entities edge, face and cell [30,31], along with the nodes from the primal grid. These partition general polyhedra into tetrahedra, where every tetrahedron is made up of a node, edge, face and cell. Edge is the midpoint between two nodes, face is the average position of the nodes describing a polyhedral face, and cell is the average position of all nodes of the polyhedron, see Fig. 1(a). A 2D variant of this structure including the narrow-band is seen in Fig. 1(b). We will refer to this tetrahedral grid as the ‘subgrid’. n

3. Interpolation The present work requires cubic Hermite interpolants for 1-, 2- and 3-simplices, i.e., line segments, triangles and tetrahedra. While a derivation of these are given in Appendix B, we introduce the concepts and methods required to evaluate the interpolants here. It is advantageous to work with barycentric coordinates, which for tetrahedra are defined by

p = T λ,

 

x , p= 1



x1 T = 1

x2 1

x3 1



x4 , 1

⎡

⎤ λ1 ⎢λ ⎥ λ = ⎣ 2 ⎦, λ3 λ4

x=

x1 x2 x3

,

(3)

where xi is the position of vertex i, x is an arbitrary point, and λ are the corresponding barycentric coordinates. Conversely,

λ = T −1 p .

(4)

50

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

Similar expressions hold mutatis mutandis for the barycentric coordinates of triangles and line segments. Barycentric coordinates are often referred to as volume (tetrahedron) or area (triangle) coordinates. For a point x, the physical interpretation of the barycentric coordinate λi pertaining to vertex i is the ratio of the volume of the tetrahedron formed by x and all vertices except i to the volume of the whole tetrahedron, assuming that x is located inside the tetrahedron. Thus, the four barycentric coordinates describing any point x sum to one, which is also seen from Eq. (3). While we will derive a cubic interpolant, the barycentric coordinates are in fact exactly the vertex weights for linear interpolation. A function f (λ) given in barycentric coordinates can be differentiated in Cartesian coordinates using the chain rule. Using Einstein’s summation convention,

∂ f (λ) −1 ∂ f (λ) ∂ f (λ) ∂λm = = T , ∂ xi ∂λm ∂ xi ∂λm mi ∂ f (λ) −1 −1 ∂ f (λ) = T T . ∂ xi ∂ x j ∂λm ∂λn mi nj

(5)

In the following, the summation convention is only invoked when stated explicitly. For more information on the topics of this section, the reader should consult e.g. [32]. 3.1. Cubic Hermite interpolation on simplices We express the cubic Hermite interpolants for simplices in Bernstein–Bézier form (B-form). This simplifies several procedures associated with the interpolant, such as differentiation and geometric interpretation. Moreover, all of the Bernstein basis polynomials are linearly independent. The Bernstein basis polynomials and corresponding coefficients employ a multiindex notation, where each index is associated with a vertex of the simplex. The index-combinations that are not given explicitly below are found by permutation of the multi-index, which merely corresponds to rearranging the vertex numbering. For example, for the tetrahedral interpolant, c 2001 may be found by swapping the second and fourth index on both sides of the equation for c 2100 . For notational convenience, we have also assigned a multi-index to the following vertex data: position x, function value f and gradient g. These should be interpreted as barycentric coordinates – e.g. x010 is the position of vertex 2 in a triangle. The 1- 2- and 3-simplex interpolants are given below. 3.1.1. 1-simplex interpolant (line segment)



P (λ) =

c i j B 3i j ,

i + j =3

B 3i j =

3!

j

λ1i λ2 ,

i! j! c 30 = f 10 ,

c 21 = f 10 +

1 3

(x01 − x10 ) · g 10 .

(6)

3.1.2. 2-simplex interpolant (triangle)



P (λ) =

c i jk B 3i jk ,

i + j +k=3

B 3i jk =

3!

j

i ! j !k! c 300 = f 100 ,

λ1i λ2 λk3 , 1

c 210 = f 100 + c 111 =

1 4

3

(x010 − x100 ) · g 100 , 1

(c 210 + c 201 + c 120 + c 102 + c 021 + c 012 ) − (c 300 + c 030 + c 003 ). 6

3.1.3. 3-simplex interpolant (tetrahedron)

P (λ) =



c i jkl B 3i jkl ,

i + j +k+l=3

B 3i jkl =

3!

i ! j !k!l! c 3000 = f 1000 ,

j

λ1i λ2 λk3 λl4 ,

(7)

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

c 2100 = f 1000 + c 1110 =

1 4

1 3

51

(x0100 − x1000 ) · g 1000 , 1

(c 2100 + c 2010 + c 1200 + c 1020 + c 0210 + c 0120 ) − (c 3000 + c 0300 + c 0030 ). 6

(8)

3.2. Gradient limiting To suppress overshoots, a gradient limiter is imposed on the interpolants given by Eqs. (6)–(8). The importance of gradient limiting is further discussed and demonstrated in Section 6. It changes the gradients locally on a ‘per simplex’ basis, and thus breaks the nodal C 1 continuity. The gradient limiter is applied edge-wise to all edges in the simplex. We will use the multi-index notation of the 1-simplex, noting that the edges of 2- and 3-simplices may be identified and reduced to two indices by the following procedure: For an m-simplex, remove m − 1 of the entries in the multi-index, e.g. j and l in the multi-index i jkl. Then, exactly four of the Bézier-coefficients will have multi-indices 30, 21, 12 and 03. These four determine the interpolant along an edge of the simplex; if j and l are removed as in the above example, the edge is between vertex 1 and 3. The other edges are obtained by all permutations of the two removed indices. For notational convenience, we define the following differences:

D 1 = 3(c 21 − c 30 ), D 2 = 3(c 03 − c 12 ), 3 D 3 = (c 12 − c 30 ), 2 3 D 4 = (c 03 − c 21 ), 2 D 5 = c 03 − c 30 .

(9)

The gradient limiting procedure is given in Algorithm 1. The Bézier coefficients c 21 and c 12 affect the end-point gradients, but not the function value. Altering these two associated with one edge does not affect the gradient along any of the other edges, where the associated basis polynomials equal zero. Algorithm 1 may thus be applied sequentially. The Bézier coefficients with three ones in the multi-index should be computed after limiting the gradient, as they are given in terms of potentially gradient-limited Bézier coefficients. Algorithm 1 Gradient limiting. if ( D 1 − D 3 )( D 1 − D 5 ) < 0 then c 12 ← 2c 21 − c 30 else if ( D 2 − D 4 )( D 2 − D 5 ) < 0 then c 21 ← 2c 12 − c 03 else if ( D 1 − D 5 )( D 2 − D 5 )  0 then c 21 ← 23 c 30 + 13 c 03 c 12 ← end if

2 c 3 03

+ 13 c 30

While the form of the gradient limiter is partially heuristic, it is based on considerations that arise from the study of one-dimensional distance functions (see Sections 5 and B.4). When these are discretely sampled and approximated by cubic Hermite interpolation, overshoots are produced if the function is sampled near a gradient discontinuity. In level set methods, these occur close to the interface for poorly resolved geometrical features, and in order not to produce additional mass, overshoots should be suppressed. An illustrating example is given in Section 6.2. 4. Semi-Lagrangian transport A semi-Lagrangian formulation is used to advect the level set scalar and gradient. While a fully Lagrangian method stores the data in arbitrary points following the velocity field, a semi-Lagrangian method stores the data in stationary computational nodes. The advection is carried out by backtracing pathlines. We transport both a scalar function and its gradient. The tetrahedral cubic Hermite interpolant from the previous section is used to approximate the scalar and gradient field between nodes. The passive scalar advection equation reads

DΦ Dt D Dt

= 0,

=

∂ + u · ∇. ∂t

(10)

D In the above equation, Dt denotes the material derivative, i.e., the time derivative following a material element, and u denotes the advective velocity field. Taking the gradient of Eq. (10), one obtains the advection equation for the gradient:

52

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

DΨ

= −∇ u · Ψ ,

Dt

Ψ = ∇Φ.

(11)

Eqs. (10) and (11) describe time derivatives along the pathlines given by

Dx Dt

  = u x(t ), t .

(12)

At this point, it is useful to anchor the material point x some place in space and time. Following [15], we thus replace x(t ) by X (x, t , τ ), where x denotes the point on the pathline X corresponding to time t, and τ is a free time variable providing all points on the pathline. Thus,

dX dτ

  = u X (x, t , τ ), τ ,

X (x, t , t ) = x.

(13)

From Eq. (13), one can find the point of departure. That is, the material point xd whose position, when advected by u, coincides with the arbitrarily chosen point x after a time step:

xd = X (x, t + t , t ).

(14)

Advancing Eq. (10) one time step t amounts to evaluating

Φ(x, t + t ) = Φ(xd , t ).

(15)

The time-update for Φ in a grid node located at x g is thus found by choosing X (x g , t + t , t + t ) = x g and integrating Eq. (13) backwards numerically, from τ = t + t to τ = t, followed by an evaluation of Eq. (15). In Eq. (11), there is a source term that changes along the pathline. While Eq. (11) can be solved directly to transport the gradient, one can also take the gradient of the discretized scalar transport equation (15) to obtain a simpler update rule for the gradient. Noting that the point of departure is a function of x, one can use the chain rule to express the gradient of the scalar function at the next time step in terms of the gradient at the point of departure:

Ψ (x, t + t ) = ∇ xd · Ψ (xd , t ).

(16)

This is referred to in [15] as a superconsistent scheme, providing better coherence between the advection of Ψ and Φ , and it has the additional advantage that Ψ only needs to be evaluated once per time step, regardless of the chosen time integration scheme. The term ∇ xd however, depends on the time integration scheme. Like [15], we use a third order Total Variation Diminishing Runge–Kutta scheme (TVD RK3) [33] for the solution of Eq. (13):

x1 = x − t u (x, t + t ), x2 =

3

x+

1

x1 −

1

t u (x1 , t ),   1 2 2 1 xd ≈ x + x2 − t u x2 , t + t . 4

3

4

3

4

3

(17)

2

Using the TVD RK3 scheme, ∇ xd in Eq. (16) is found as

∇ x1 = I − t ∇ u (x, t + t ), ∇ x2 =

3

3

I+

1

1

∇ x1 − t ∇ x1 · ∇ u (x1 , t ), 4   1 2 2 1 ∇ xd ≈ I + ∇ x2 − t ∇ x2 · ∇ u x2 , t + t . 4

4

3

3

2

(18)

∇ xd can be computed in parallel with xd in Eq. (17). Unless the velocity field is given analytically, the velocity and velocity gradient must be interpolated in space and time from the nodal values. The procedure for advancing Φ and Ψ one time step in a grid node is summarized in Algorithm 2.

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

53

Algorithm 2 Advection – updating grid node g from time t to t + t.

x ← xg . Solve Eqs. (17)–(18). Find the tetrahedron T in which xd is located, and number its vertices from 1 to 4. Find the barycentric coordinates λ of xd with respect to T from Eqs. (3)–(4). Evaluate Eqs. (15)–(16). Approximate Φ(xd , t ) by the interpolant P (λ) from Eq. (8), subject to Algorithm 1. Approximate ∇Φ(xd , t ) by Eq. (5) applied to P (λ).

5. Reinitialization of the level set function When advecting an initially smooth level set function, it becomes distorted, and the regions containing the interface may be smeared out or strongly compressed. To maintain a numerically well-defined interface, the initial smoothness thus needs to be restored from time to time. The process of restoring the initial smoothness is usually referred to as reinitialization in the level set literature. It is of great importance that the reinitialization minimally disturbs the zero-level set, lest the interface position is altered. The most common approach is to let the level set function away from the interface approximate a signed distance function, i.e., a function whose absolute value is the distance to the closest point on the interface (zero-level set), and whose sign changes as one crosses the interface. For this purpose, there exist several methods. The most common is the solution of an eikonal equation, which describes the travel time of an advancing front emanating from the interface:

∇Φ = 1, Φ(xΓ ) = 0,

∀xΓ ∈ Γ,

(19)

where Γ is the set of points belonging to the interface. To obtain a single-valued solution of Eq. (19) corresponding to the first arrival time, a small diffusive term ∇ 2 Φ can be added to the right hand side [34], which yields the ‘viscosity solution’ [35] as approaches zero. To account for the sign change across the interface,









sgn Φ(x) = sgn Φ0 (x) ,

(20)

where sgn is the signum function, and Φ0 is the initial level set function from which Γ is found. Assuming a uniform propagation of unit velocity, the travel time of the advancing front to reach any point in the domain gives the distance to the closest point on the interface. While the eikonal equation may be treated as a time dependent PDE and solved iteratively from an initial state [36], it may also be solved in a single pass due to the upwind-nature of the front propagation problem. This requires a strict adherence to the upwind causality, which can be enforced through a priority queue of the grid nodes right outside the front. Such a method is called a Fast Marching Method (FMM), as it ‘marches’ the solution outwards from the interface. It is similar in nature to Dijkstra’s algorithm [37], but whereas the latter is restricted to a path following the nodal connectivity, an FMM approximates the Euclidean distance. The remaining part of this section is structured as follows. The methodology of the FMM is given in Section 5.1, the local travel time approximation is described in Section 5.2, the initialization of the FMM is described in Section 5.3, a discussion of the expected accuracy is given in Section 5.4, and a brief explanation of how the algorithm can be parallelized is given in Section 5.5. More details on the method and its relation to the eikonal equation are given in Appendix A. 5.1. The Fast Marching Method As the Fast Marching Method sequentially ‘marches’ the solution outwards from the interface, a book-keeping system is needed to keep track of which nodes are already processed and which nodes are prospective candidates for update, i.e., which nodes lie close to the current front position. Grid nodes are classified into three categories: accepted : close : far :

Nodes that lie inside the front, and have received their final travel time/distance value. Nodes that lie outside, but adjacent to the front. Nodes that lie outside the front, and are not in the close set.

In the current setting, ‘adjacent’ means immediate neighbors on the grid. Note that the names of the three sets vary in the literature. To ensure the data transfer is a strict upwinding process, the close set is sorted according to tentative travel times computed from the accepted set. The sorting places the nodes in a priority queue, and at every step of the algorithm, the close node with the smallest tentative value (i.e., the node that is first encountered by the moving front), is removed from the close set and added to the accepted set. All of its neighbors that were in the far set are then transfered to the close set. Tentative values are computed for these, and if some of the neighbors were in the original close set, their values are recomputed to take into account the newly available information. The process is repeated until all of the nodes are in the accepted set. The sorting process associated with the close nodes is a potential bottleneck. But since few changes are made to the order in a single step of the algorithm, an effective binary tree structure known as a ‘heap’ may be used. For detailed descriptions of the heap structure and associated operations, the reader is referred to [38]. For a comprehensive review of Fast Marching Methods, see [34].

54

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

Fig. 2. Accepted and close nodes.

5.2. Local travel time approximation In the following discussion, it is important to keep in mind the equivalence between travel time and distance due to the unit propagation velocity, as the two will be used interchangeably. The level set function is in our method always discretized on a tetrahedral grid. A node in the close set may have either 1, 2 or 3 local neighbors in an adjacent tetrahedron that belong to the accepted set. By adjacent, we mean a tetrahedron in which this node is a vertex, and we will refer to it as the ‘close node’. The accepted nodes in an adjacent tetrahedron will form either a point, a line segment or a triangle (0-, 1- or 2-simplex respectively), and in the cases of 2 or 3 neighbors, we assume that the accepted solution may be interpolated on the simplex formed by the accepted neighbors. We will refer to such a simplex as an ‘accepted simplex’, and its convex hull is denoted by T . An example is given in Fig. 2. Here, the close node C 1 has two accepted 2-simplices; ( A1, A2, A3) and ( A2, A3, A4). The close node C 2 has only one accepted 1-simplex ( A3, A4). We use Huygens’ principle [39] as formulated in [40] to determine the value of the distance function at the close node: ‘. . . the wave front is the envelope of secondary waves whose centers are themselves on a previous wave front.’ With a constant propagation velocity, the secondary wave fronts are spherical. Huygens’ principle has difficulties explaining the implied presence of backward-traveling waves, and they are thus neglected. However, our sole interest in the first arrival time makes the use of Huygens’ principle unproblematic, and the idea of its use with level sets is not new. The arrival time f at the close node for the spherical front that is emitted from a point x in the accepted simplex, is the time at which the spherical front is emitted plus the distance from the emitter, as the front propagates with unit velocity. Using an interpolant P to approximate the distance function in T , the arrival time can be estimated as

f (x) ≈ P (x) + D (x),

(21) c

where D (x) is the Euclidean distance from x to the close node’s position x . The envelope of such fronts describes a section of the main front emanating from the accepted simplex. To obtain the first arrival time at the close node of the fronts emanating from T , one must find the point x in T that minimizes f . To consider this point valid in updating the close node, it must satisfy two criteria: 1. x must lie in T , i.e., all of its barycentric coordinates with respect to the accepted simplex must be nonnegative. 2. f (x) must be greater than the distance function in all of the accepted nodes from which it is computed. The first criterion is to avoid extrapolation into neighboring tetrahedra, where a local interpolant should be used instead. The second criterion ensures that the close node is updated from upwind values only. To obtain the first arrival time at xc of the fronts emanating from all adjacent accepted simplices, we solve

F = min f A (x), A , x∈ T A

(22)

where T A is the convex hull of an adjacent accepted simplex A, and f A includes the appropriate expressions for P and D on A. Assuming that only one point in any T A admits the minimum value F , we define xe to be this point, i.e., the point that emits the spherical wavelet that first reaches the close node:

xe = arg min f A (x). A , x∈ T A

In practice, if by chance several points provide F , the implemented code picks one based on the node indexing. represents a constrained minimization problem. Using an unconstrained minimization algorithm, one must search minima in A, and all of its faces. A face is in this respect any simplex of lower dimension that can be formed from of the vertices in A. If A is a triangle, the search for local minima must thus include the edges of the triangle

(23) Eq. (22) for local a subset and the

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

55

Fig. 3. Variables and local coordinate system (see Appendix A) used to determine the travel time in a tetrahedron. A spherical wavelet emanates from x, and reaches the close node in D time.

vertices, and if A is a line segment, the search must also include the vertices (end-points). It is important to note that since f is not necessarily a convex function, the minimization algorithm may fail to converge to the true minimum. A good initial guess is thus important, see Appendix A. To facilitate the use of cubic Hermite interpolation, we also store the distance function’s gradient. When the point of emission xe providing F is found, tentative values of Φ c and the gradient Ψ c of the level set function in the close node are updated according to

Φc = F , Ψc =

xc − xe

xc − xe 

,

(24)

where xc is the position vector of the close node. These are not terminal values unless Φ c is the smallest value of Φ in the close set when all close nodes have been updated in this manner. In that case, the close node is transfered to the accepted set. The update rule for Ψ in Eq. (24) is a consequence of the properties of the distance function; its gradient has unit magnitude, and the gradient’s direction is normal to the propagating front. Our approach is illustrated in Fig. 3. Similar update rules are found in [23,24]. In [41], Fermat’s principle was applied to arrive at the same type of update rule, and it was shown that with a linear interpolant on Cartesian grids, it reduces to the update rule obtained by FD discretization of the eikonal equation. The work of [41] was extended to general unstructured tetrahedral grids in [42]. The present work may be seen as an extension of [42] to cubic Hermite interpolation stencils. Note the equivalence between Fermat’s principle and Huygens’ principle [43], and that the FMM itself is a consequence of Huygens’ principle [44]. In the present work, we advance both the distance function and its gradient for use in the cubic Hermite interpolant. While improving the accuracy, this necessitates the use of a robust optimization algorithm; f may in fact have several local minima in high-curvature regions and near colliding fronts. As both the gradient and Hessian of P and D are easily computed, we use Newton’s method with a trust-region approach. We use the algorithm library known as TOMS611 [45]1 as a ‘black-box’ solver, employing the HUMSL subroutine described in [45]. The iterative method needs an initial guess for x which is readily provided by using linear rather than cubic Hermite interpolation on the accepted simplex. Linear interpolation provides a unique solution unless the linear gradient has a magnitude larger than unity. Otherwise, the projection of the close node down on the accepted simplex is used as an initial guess. See Appendix A for details. Where two advancing fronts meet, the gradient of the distance function is discontinuous. In high-curvature regions and thin films, this happens near the interface. Clearly, a polynomial interpolant is not well suited here, and the order of accuracy is reduced, as noted in [15]. As demonstrated in Section 6.2, the FMM is especially prone to overshoots in these areas, and the interpolants given in Eqs. (6)–(7) are thus subject to the gradient limiting of Algorithm 1. Finally, we note that the requirement of compact stencils has one shortcoming in the FMM context – there might not exist a sufficient number of upwind nodes for a proper 3D reconstruction of the advancing front. The problem is illustrated in 2D in [46], where it is shown that with a linear front, a strictly acute triangulation is required to ensure a sufficient number of upwind nodes with compact stencils. Alternatively, a remeshing method is proposed in [46], to ensure acute angles. In [42], the upwind requirement is instead relaxed, noting that this might be an acceptable approach if the front is sufficiently smooth. In the present work, we reconstruct the front from one or two nodes if a proper reconstruction cannot be done from three accepted nodes. We consider this an acceptable approach, even though O (h) and O (1) errors may be locally introduced in Φ and Ψ respectively, where h is the grid spacing. As the test cases in Section 6 demonstrate, accurate solutions and good convergence in the L 1 error norm are still obtained.

1

Source code available on Netlib: http://www.netlib.org/toms/611.

56

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

5.3. Initialization An initial accepted set is required to start the FMM. Rather than labeling a set of nodes as accepted – as in e.g. [47] – we use edge intersections of the zero-isosurface. The edge intersections replace neighbor nodes in the connectivity graph. As an example, if an edge intersection is found between node 1 and 2, node 2 is replaced by the edge intersection in node 1’s neighbor list, and vice versa. The edge intersections constitute the initial accepted set, and carry information on position, gradient (interface normal) and level set value; the latter is obviously always zero. This requires keeping a list of interface intersections and corresponding interface normals. In the FMM routine, every time the code attempts to access neighbor information on a node that is located on the other side of the interface, it should be redirected to the edge intersection list. The edge intersections partition the domain into two parts; one positive and one negative, between which information cannot be exchanged. Hence, to avoid the clutter of locally applying the correct sign in the numerical procedures, the positive and negative sides are flagged initially, and the initial interface normals on the negative side are reversed. The positive and negative sides are then treated in two separate passes using the same ‘positive’ discretization, before the signs of function value and derivatives on the negative side are reverted. To obtain the edge intersections, the cubic Hermite interpolant of Eq. (6) is used. This results in a cubic expression whose real roots describe the positions of the edge intersections. An initial gradient is also required. The gradient component along the edge is found by differentiating the interpolant, while the edge-orthogonal components are found by linear interpolation to the point of intersection. The two are then added, and the resulting gradient vector is normalized to account for the distance function property. The gradient limiting of Algorithm 1 is applied to the cubic Hermite interpolant before solution of the cubic equation. Initialization of an edge between nodes i and j is algorithmically summarized in Algorithm 3. Algorithm 3 Initialization algorithm. For an edge between nodes i and j, label node i as vertex 1 and node j as vertex 2. Obtain the Bézier coefficients from Eq. (6). Subject the Bézier coefficients to Algorithm 1. a ← 3c 21 − c 30 − 3c 12 + c 03 b ← 3c 30 + 3c 12 − 6c 21 c ← −3c 30 + 3c 21 d ← c 30 T ← {t ∈ R: 0  t  1 ∧ P (t ) = at 3 + bt 2 + ct + d = 0} if T = ∅ then t 1 ← inf( T ) t 2 ← sup( T ) for k = 1, 2 do xkint ← tk x j + (1 − tk )xi x −x

v 1 ← x j−x i2 P  (tk ) j i v 2 ← tk Ψ j + (1 − tk )Ψ i v3 ← v1 + v2 −

Ψ kint ←

((x j −xi )· v 2 )(x j −xi ) x j −xi 2

v3

v3 

end for int Replace the position and gradient in node j as a neighbor of i by xint 1 and Ψ 1 .

int Replace the position and gradient in node i as a neighbor of j by xint 2 and Ψ 2 . else return end if return

In Algorithm 3, t is λ2 , obtained by substituting λ1 = 1 − λ2 . t = 0 thus gives the position of node 1, and t = 1 gives the position of node 2. As can be seen, if several zero-contour crossings are found along an edge, the ones closest to each node on either side are used to replace the neighbor connectivity. 5.4. Expected accuracy Assuming a smooth distance function, a sufficient number of upwind neighbors and the use of the cubic Hermite interpolant of Eq. (7), we expect the local approximation of Φ according to Eq. (22) to introduce O (h3 ) errors due to the truncation error in Eq. (B.4). This means that the estimated point of emission xe will be offset from the true point of emission by an O (h3 ) distance, where h is the grid spacing. From this, we can determine the expected accuracy of Ψ , found by Eq. (24). As h decreases towards zero, the computed gradient will be rotated an angle of O (h2 ) as compared to the true gradient. This leaves the gradient with an O (h2 ) error. The FMM can overall not be more accurate than the initial accepted set, as errors in this set will be propagated to the entire domain. As the univariate cubic Hermite interpolant is used to determine the initial zero-level set, the errors introduced in Φ are O (h4 ). The linear interpolation used to find the edgeorthogonal gradient components introduces an O (h2 ) error in Ψ . We thus expect, under optimal conditions, convergence rates of third order in Φ , and second order in Ψ in a band of constant thickness with respect to the number of cells around

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

57

the interface. This may however be reduced by non-smooth geometry and grids with highly obtuse angles, as noted in Section 5.2. If the FMM is advanced in the whole domain, a reduction of one in the order of convergence can be expected. This is because the number of cells away from the interface increases by h−1 , which is multiplied by the local discretization error to take error accumulation into account. 5.5. Parallelization While the FMM is in essence a serial algorithm, it is clear that different parts of the front may be advanced in parallel, as two parts of the front located some distance apart are not mutually dependent until they eventually merge. Spatially partitioning the advancement of the FMM into several processes, the processes may run in parallel, while the FMM is executed serially on each process. Upwind conflicts may thus arise. When process A receives a value of Φ from another process B at a shared node, this value must be compared to the largest accepted value existing on process A. If the received value is smaller, some of the accepted values on process A may not have been computed based on the correct causality, and the solution must be stepwise reverted on process A until the value received from B is larger than the largest accepted value on A. A parallel method based on the approach of [48] is implemented with the present discretization. While we will not go into further detail on parallel implementation and performance, we note that the only amendment required to the method of [48] is the transfer of gradients at shared nodes. 6. Numerical tests In this section we test the method’s performance. As error indicators in test cases where the primary interest is interface motion, we use volume and shape preservation. In cases where the aim is to assess the accuracy everywhere in the computational domain we use the L 1 and L ∞ error norms of Φ and Ψ . The total positive volume is approximated by applying the Heaviside function to a piecewise linear reconstruction of Φ in each tetrahedron. The Hermite interpolant is first used to assign values of Φ to a grid that is roughly four times finer, applying the partitioning method described in Section 2 recursively. The volume of the positive phase is then found by





V+ =



H Φ(x) dx,

(25)

Ω

where the integral is approximated by a linear, geometric reconstruction of the interface. As a measure of the average error in interface position, we use the error measure suggested in [49]:



1

EΓ =

A

  ex     H Φ (x) − H Φ comp (x)  dx,

(26)

Ω

where Φ ex denotes the exact solution and Φ comp denotes the computed solution. The integral in Eq. (26) over one tetracomp ex hedron is approximated by | V + − V + |, where V + is computed as in Eq. (25), but the integration is constrained to the tetrahedron. This simplification under-predicts the error when the computed and exact interface intersect. Sensitivity studies have shown that this effect is insignificant. We define the error in Φ and Ψ in node i as









i eΦ = Φ ex (xi ) − Φ comp (xi ),

i eΨ = Ψ ex (xi ) − Ψ comp (xi ).

(27)

The L 1 and L ∞ error norms are evaluated according to

L1 =

1  N

ei ,

i

L ∞ = max e i , i

(28)

where N is the number of nodes under consideration. In the following tests, we use the exact analytic expressions for the velocity and velocity derivatives to solve the transport equations for the level set scalar and gradient. 6.1. Deformation field test – convergence To assess the global convergence rate of the advection scheme, we use the velocity field suggested in [50,51] with the smooth scalar field used in [15]. As the scalar field is smooth, there should be no need for gradient limiting. However, as the smooth function is distorted by the velocity field, the discretized field might not be perceived as smooth everywhere, and

58

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

Table 1 Results from the deformation field test, with gradient limiting. h

L 1 error in Φ

p

L 1 error in Ψ

p

L ∞ error in Φ

p

L ∞ error in Ψ

p

1/16 1/32 1/64 1/128

2.09 × 10−4 2.87 × 10−5 3.84 × 10−6 4.92 × 10−7

– 2.9 2.9 3.0

6.05 × 10−3 1.06 × 10−3 1.57 × 10−4 2.35 × 10−5

– 2.5 2.8 2.7

1.36 × 10−3 2.05 × 10−4 3.50 × 10−5 5.81 × 10−6

– 2.7 2.6 2.6

1.06 × 10−1 2.79 × 10−2 1.28 × 10−2 9.05 × 10−4

– 1.9 1.1 3.8

Table 2 Results from the deformation field test, without gradient limiting. h

L 1 error in Φ

p

L 1 error in Ψ

p

L ∞ error in Φ

p

L ∞ error in Ψ

p

1/16 1/32 1/64 1/128

1.76 × 10−4 2.64 × 10−5 3.62 × 10−6 4.67 × 10−7

– 2.7 2.9 3.0

4.01 × 10−3 7.38 × 10−4 1.13 × 10−4 1.62 × 10−5

– 2.4 2.7 2.8

9.49 × 10−4 1.61 × 10−4 2.31 × 10−5 3.06 × 10−6

– 2.6 2.8 2.9

1.62 × 10−2 4.54 × 10−3 9.85 × 10−4 2.97 × 10−4

– 1.8 2.2 1.7

Fig. 4. Equally spaced isocontours of Φ at t = 1.

gradient limiting might be triggered by Algorithm 1. We include cases both with and without application of Algorithm 1, to check its effect on the accuracy and convergence rate. The velocity field is 2D, and given by the stream function

ψ(x, y ) =

1

π

 sin2 (π x) sin2 (π y ) cos

πt T

 ,

(29)

with T = 2. The resulting velocity field is constrained to a [0, 1] × [0, 1] box, with vanishing velocity on the boundaries. The advected scalar should return to its initial condition at t = 2, which is where we measure the error. The initial scalar field is given by

     Φ(x, y ) = exp − (x − x0 )2 + ( y − y 0 )2 − exp −r02 , x0 = 0.5, y 0 = 0.75, r0 = 0.15.

(30)

Fig. 4 shows the scalar field in its most distorted state. We use a 3D discretization for a 2D problem, and the domain is thus given unit length over a single cell in the z-direction, with periodic boundary conditions. A constant time step of t = 0.002 is used with Cartesian primal grids of spacings h = 1/16, h = 1/32, h = 1/64 and h = 1/128. The resulting tetrahedral subgrid is twice as fine. The results with and without gradient limiting are summarized in Tables 1 and 2 respectively, where p denotes the estimated order of convergence. Without gradient limiting, the order of convergence is around three for Φ and two for Ψ in the L ∞ error norm. This is as expected – the triangular and tetrahedral cubic Hermite interpolants’ accuracy are limited by the approximation of the ‘bubble functions’ by Eqs. (B.4)–(B.7). While Φ shows the same convergence in the L 1 error norm, Ψ converges slightly faster than expected. While we have not found an obvious reason, this might be caused by local error cancellation effects as the velocity field is reversed. With the gradient limiting algorithm applied, the L 1 error norms converge at approximately the same rates as without gradient limiting. But the L ∞ error norm of Ψ shows a less predictable behavior, as expected.

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

59

Fig. 5. Top row: Interface at t = 0 (black) and t = 8 (red). Bottom row: Interface at t = 4. Left to right: Case 1, 2 and 3. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

6.2. Deformation field test – behavior with under-resolved interface To test how the proposed method behaves when the interface geometry is insufficiently resolved, the deformation field test case was is run with T = 8. This stretches the interface much more, and using a 64 × 64 subgrid, the results are directly comparable to [15]. Three tests are run: one without reinitialization or gradient limiting (Case 1), one with reinitialization every time step and gradient limiting (Case 2), and one with reinitialization every time step, but without gradient limiting (Case 3). The results are visualized2 in Fig. 5. With a 64 × 64 subgrid, it is clear that the thinnest parts of the spiral are too thin to be properly resolved on the grid. If the FMM cannot find an interface crossing along an edge, the nodes of this edge will have the level set function overwritten by characteristics coming other parts of the interface. Thus, subgrid features are filtered out with the FMM, while they are retained (although they may be too small to produce interface intersections) without the FMM. When the velocity field is reversed, subgrid features which do not contribute to the interface in terms of grid intersections may grow back to a size where they contribute to the interface. With an FMM, they are permanently lost. As can be seen from Fig. 5, the case without reinitialization gives the best results. Compared to [15], the results are similar. When reinitialization is applied, the results are more similar to those produced by the fifth order WENO scheme implemented in [15], which includes reinitialization. Here, the mass loss is substantial, which is a characteristic of nonconservative level set methods when one phase is stretched into features that are too small to be resolved on the grid. A physical example of this is atomization, for which conservative schemes should be used [11]. Finally, reinitialization without gradient limiting produces an excessive amount of mass, giving practically useless results. The effect is especially pronounced in thin regions or near high interface curvature (see also Section 6.5). The deformation field test case thus provides a proper stress test of the gradient limiting algorithm. 6.3. Computing the signed distance function of a spherical interface To assess the accuracy of the proposed FMM, we use it to compute the signed distance and its gradient field from a spherical interface at different grid resolutions. This problem has a very simple analytical solution:

Φ = R − x0 − x, x0 − x Ψ= , x0 − x

(31)

where R is the radius of the sphere and x0 is its center point. We choose R = 0.15 and x0 = (0.5, 0.5, 0.5) in a [0, 1] × [0, 1] × [0, 1] domain with a Cartesian primary grid. Three sets of convergence tests are run: 2

The apparent jaggedness is an artifact from the linear interpolation used in the visualization software.

60

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

Table 3 Convergence results from set 1. h

L 1 error in Φ

p

L 1 error in Ψ

p

L ∞ error in Φ

p

L ∞ error in Ψ

p

1/8 1/16 1/32 1/64 1/128

8.70 × 10−5 1.01 × 10−5 1.37 × 10−6 4.22 × 10−7 1.26 × 10−7

– 3.1 2.9 1.7 1.7

2.49 × 10−3 4.47 × 10−4 1.04 × 10−4 2.80 × 10−5 8.00 × 10−6

– 2.5 2.1 1.9 1.8

3.01 × 10−4 6.87 × 10−5 3.29 × 10−5 1.62 × 10−5 7.99 × 10−6

– 2.1 1.0 1.0 1.0

9.21 × 10−3 8.29 × 10−3 1.03 × 10−2 1.11 × 10−2 1.07 × 10−2

– 0.2 − 0.3 − 0.1 0.1

Table 4 Convergence results from set 2. h

L 1 error in Φ

p

L 1 error in Ψ

p

L ∞ error in Φ

p

L ∞ error in Ψ

p

1/8 1/16 1/32 1/64 1/128

9.79 × 10−3 5.02 × 10−3 2.50 × 10−3 1.26 × 10−3 6.30 × 10−4

– 1.0 1.0 1.0 1.0

4.48 × 10−2 2.58 × 10−2 1.37 × 10−2 7.12 × 10−3 3.65 × 10−3

– 0.8 0.9 0.9 1.0

1.72 × 10−2 9.22 × 10−3 5.62 × 10−3 3.48 × 10−3 2.11 × 10−3

– 0.9 0.7 0.7 0.7

4.63 × 10−1 4.72 × 10−1 4.76 × 10−1 4.79 × 10−1 4.82 × 10−1

– 0.0 0.0 0.0 0.0

Table 5 Convergence results from set 3. h

L 1 error in Φ

p

L 1 error in Ψ

p

L ∞ error in Φ

p

L ∞ error in Ψ

p

1/8 1/16 1/32 1/64 1/128

1.13 × 10−5 8.61 × 10−7 9.90 × 10−8 1.23 × 10−8 1.56 × 10−9

– 3.7 3.1 3.0 3.0

6.95 × 10−4 1.46 × 10−4 3.96 × 10−5 1.07 × 10−5 2.81 × 10−6

– 2.3 1.9 1.9 1.9

5.11 × 10−4 4.57 × 10−5 2.28 × 10−5 1.14 × 10−5 5.71 × 10−6

– 3.5 1.0 1.0 1.0

3.33 × 10−2 8.60 × 10−3 8.60 × 10−3 8.60 × 10−3 8.60 × 10−3

– 2.0 0.0 0.0 0.0

Table 6 Convergence results from the FMM initialization procedure. h

L 1 error in Φ

p

L 1 error in Ψ

p

L ∞ error in Φ

p

L ∞ error in Ψ

p

1/8 1/16 1/32 1/64 1/128

4.87 × 10−5 2.34 × 10−6 1.67 × 10−7 9.43 × 10−9 5.72 × 10−10

– 4.4 3.8 4.1 4.0

6.82 × 10−3 2.17 × 10−3 5.31 × 10−4 1.22 × 10−4 2.95 × 10−5

– 1.7 2.0 2.1 2.0

1.80 × 10−4 1.24 × 10−5 9.67 × 10−7 6.23 × 10−8 3.97 × 10−9

– 3.9 3.7 4.0 4.0

2.48 × 10−2 6.20 × 10−3 2.19 × 10−3 5.62 × 10−4 1.42 × 10−4

– 2.0 1.5 2.0 2.0

Set 1 To initialize the FMM, the analytical solution is prescribed at every grid point. This only serves to obtain the initial accepted set, which is found by Algorithm 3. The update procedure proceeds as described in Section 5. Cubic Hermite interpolation is used to approximate the solution on accepted simplices. Set 2 As set 1, but with linear interpolation approximating the solution on accepted simplices. Set 3 The whole FMM procedure is discarded to isolate the effect of the local travel time approximation. Cubic Hermite interpolation is used to approximate the solution on accepted simplices, which have vertex values set to the analytical solution. Thus, no initialization is required. While we expect third order convergence in Φ in the L 1 error norm in the third set of tests, we expect slower convergence in the first set of tests. This is because error accumulates as the front propagates outwards. To reach a point a distance L from the initial interface, the front must propagate across approximately L /h cells, where h is the grid spacing. As the grid is refined, more steps must be taken to reach the same point. Results from the three sets of tests are given in Tables 3–5, and the errors produced by the initialization algorithm are given in Table 6. p denotes the estimated order of convergence. In the error norms for Ψ , the center point, where Ψ is undefined, is not included. The slow convergence in the L ∞ norm of Φ may be explained by the analytical distance function. For the grid nodes closest to the center of the sphere, the third derivatives of the distance function increase as h−2 , where h is the grid spacing. In the Taylor expansion of Eq. (B.4), the leading order error term is multiplied by the third derivative, yielding first order convergence of Φ in the nodes closest to the center of the sphere. The L ∞ norm of Ψ does not converge, as expected since Ψ converges one order slower than Φ . The major difference between the approximately constant L ∞ norms

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

61

Fig. 6. Isosurfaces of Φ with the 8 × 8 × 8 primal grid, Φ = 0.1. The zero-isosurface is shaded darker. Results with linear interpolation (left), and cubic Hermite interpolation (right).

Fig. 7. Volume and interface error development for repeated reinitialization.

in Ψ between the linear and cubic Hermite interpolants should be noted. For the case with neighbors set to the analytical solution, the L 1 norm of Φ shows third order convergence in set 3, as expected from the accuracy of the triangular cubic Hermite interpolant. The effect of error accumulation is clearly visible in Table 4. While the linear interpolant introduces O (h2 ) errors in Φ locally, the L 1 error norm shows first order convergence, due to the accumulation of errors. This effect is strong with the linear interpolant, as the distance from the interface is consistently under-predicted on the concave side, and over-predicted on the convex side. This is not the case with the triangular cubic Hermite interpolant, and error accumulation does not strongly influence the convergence rate in the L 1 norm on the coarser grids. A visual comparison between the results produced by linear and cubic Hermite interpolation is given in Fig. 6, where the above-mentioned differences may be seen. It should also be noted that contrary to locally planar fronts, where an appropriate number of upwind neighbors can always be found given the absence of obtuse grid angles, this is not necessarily the case with a locally curved front. O (h) errors in Φ may thus be introduced locally, and propagated into other parts of the domain. Due to these two effects, it is difficult to predict the convergence in Φ , and we rely on empirical investigation. The convergence of Ψ appears to be much less prone to error accumulation than Φ . The initialization of the FMM converges at the expected rate. In the tests above, the narrow-band was extended to the whole domain. This has been done because the proposed FMM is applicable as a standalone solver, and may be applied to other problems than interface advection where the distance field away from the interface may be of interest. In the context of interface advection, the main interest is to reinitialize the nodes closest to the interface, and slower convergence due to error accumulation is thus not an issue. 6.4. Repeated reinitialization As reinitialization involves a resampling of the interface location, it is clear that it may introduce additional numerical diffusion. To better see the effect of frequent reinitialization, we apply the Fast Marching Method repeatedly to a stationary interface, and measure the volume and interface error as a function of the number of applications. The interface has the shape of an octant of a sphere of radius 0.3; a length which is discretized by approximately 15 subgrid spacings. The primal grid is unstructured tetrahedral. We sequentially apply 100 reinitializations and extract the interface error and ratio of current to initial volume. The resulting development is shown in Fig. 7. The initial rapid growth in error is caused by edge smoothing. After 100 iterations, around 1% of the analytical volume is lost, and the interface position is off by an average of roughly 0.2% of the radius. The initial and final shapes are shown in Fig. 8.

62

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

Fig. 8. The initial shape (left) after 100 reinitializations (right).

Table 7 Parameters in the Zalesak’s disk test. Domain Disk center Disk radius Notch with Top of notch y-coordinate u (x, y ) v (x, y )

[0, 100] × [0, 100] (×[0, 1]) (50, 75) 15 5 85

−Ω( y − 50) Ω(x − 50) π

Ω

600

Table 8 Results from Zalesak’s disk after one revolution. Gradient limiter on/off

Reinit. on/off

Area loss [%]

EΓ

E Γ order

on

off

25 × 25 50 × 50 100 × 100

−11.30 −0.32 0.13

7.38 × 10−1 8.86 × 10−2 1.58 × 10−2

– 3.1 2.5

on

25 × 25 50 × 50 100 × 100

0.41 0.07 0.10

3.10 × 10−1 6.84 × 10−2 2.27 × 10−2

– 2.2 1.6

off on

100 × 100 100 × 100

0.23 0.04

2.63 × 10−2 3.47 × 10−2

– –

off

Resolution

6.5. Zalesak’s disk Zalesak’s disk [52] is a widely used 2D advection test, testing the ability of an advection scheme to preserve the shape and area of a notched disk in a velocity field describing solid body rotation. The shape contains both smooth curves and convex/concave corners, as well as two parallel straight lines of relatively narrow spacing on coarse grids. As our proposed method is 3D in general, we use one cell of unit length in the z-direction. The geometry of the notched disk is as given in [52], and summarized in Table 7. We use a unit time step t for all spatial resolutions, and the angular frequency Ω is set such that one full revolution takes place every 1200 steps. We run tests with three different resolutions and reinitialization every time step on/off. Additionally, we test the effect of the gradient limiter by turning it off, and running the highest resolution with and without reinitialization. To avoid contamination from the fringes of the interface band, the band spans around nine primal grid cells on both sides of the interface in the cases without reinitialization. With reinitialization every time step, the band spans around two cells on either side of the interface-containing cells. The primal grid structure is unstructured prismatic triangular, and the subgrid on which the level set solution is computed is, as always, tetrahedral and roughly twice as fine. The given grid sizes refer to the total number of primal grid spacings along x- and y-aligned boundaries, and the resolution is similar in the interior of the domain. It should be noted that as opposed to a Cartesian primal grid, the prismatic triangular grid may produce tetrahedra with obtuse angles. After one revolution, we extract E Γ , which is a measure of the average error in interface position, and the change in area (volume). Both are found as described in the beginning of this section. The results are given in Table 8. Convergence without and with reinitialization is shown in Figs. 9(a) and 9(b), where the shapes with different resolutions, after one revolution, are superimposed. Fig. 9(c) shows the effect of turning off the gradient limiter, and Fig. 9(d) shows the exact solution plotted over the primal grid. As noted in [7,53], numerical diffusion in the sharply convex and concave parts approximately cancel out in terms of area preservation, so conserved area is not necessarily a good indicator of the quality of the solution. As a ‘fair’ comparison with other works, we compare our resolution on the primal grid to twice the resolution in other works, noting that the tetrahedral subgrid on which the level set method is discretized is roughly twice as fine as the primal grid. Note that additional effective resolution is introduced by the storage of gradient vectors. Compared to the results of [7],

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

63

Fig. 9. Results from Zalesak’s disk. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

the fifth order WENO scheme is outperformed, most notably on the coarsest grid, where the whole disk has vanished after one revolution. Our results are more similar to those of the Particle Level Set method where, when comparing the cases with reinitialization, we have smaller or equal interface errors, and better area preservation. The Zalesak’s disk test case was also used in [54], on a Cartesian grid with a third order WENO method for advection. Here, reinitialization was done by an FMM, which was initialized by a constrained minimization method for nodes near the interface. While E Γ was not used as an error estimate, visualization of the disk’s shape indicates accuracy similar to that of the current method. In [54], area losses are stated to be 0.652% and 0.222% for 100 × 100 and 200 × 200 grids respectively – higher than the losses obtained with the current method on the 50 × 50 and 100 × 100 grids. While Table 8 lists the accuracy, it provides little insight into the effects of the reinitialization and gradient limiter. A visual representation of the solution is therefore an important complement. Figs. 9(a) and 9(b) reveal that reinitialization serves to counteract numerical diffusion near the notch at the coarsest resolution, while it slightly increases numerical diffusion in corners at the finest resolution. It also appears to provide more symmetry about the notch for the lower resolutions. It can be seen that smoothing near corners appear to constitute most of the interface error. Here, the gradient is discontinuous and the accuracy typically drops to first order (see Section 6.3), which explains the slower convergence on finer grids. Fig. 9(c) shows the importance of gradient limiting, especially with reinitialization, where the absence of the gradient limiter results in spurious wriggles near corners. In terms of interface accuracy, both cases without gradient limiting are less accurate than their gradient limited counterparts. 6.6. Enright test With the analytical velocity field proposed in [51], the ‘Enright test’ [7], stretches a sphere into a thin film and then back again. This has become a very popular test with 3D level set methods. Depending on the spatial resolution, the stretching may or may not produce topology change, and areas of high curvature arise. We use tetrahedral primal grids,

64

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

Table 9 Parameters in the Enright test. Domain Sphere center Sphere radius u (x, y , z, t ) v (x, y , z, t ) w (x, y , z, t )

[0, 1] × [0, 1] × [0, 1] (0.35, 0.35, 0.35) 0.15 2 sin2 (π x) sin(2π y ) sin(2π z) cos( π3t ) − sin(2π x) sin2 (π y ) sin(2π z) cos( π3t ) − sin(2π x) sin(2π y ) sin2 (π z) cos( π3t )

Table 10 Results from the Enright test. Resolution

Volume loss [%]

EΓ

E Γ order

25 × 25 × 25 50 × 50 × 50

24.68 0.73

1.47 × 10−2 2.29 × 10−3

– 2.7

Fig. 10. The Enright test on an approximately 50 × 50 × 50 primal grid. From left to right: t = 0.0, t = 0.6, t = 1.5, t = 2.4, t = 3.0.

with approximate resolutions 25 × 25 × 25 and 50 × 50 × 50. The geometrical setup and velocity field are given in Table 9. The simulation is stopped at t = 3, at which the initial sphere should ideally be recovered. We use constant time steps t = 6 × 10−4 , so that the total simulation is resolved by 5000 steps. Reinitialization is applied every 200 steps, which is more realistic for practical use than every time step, as in the Zalesak’s disk test case. The bandwidth is approximately 4 cells on either side of the interface. Results are given in Table 10. At the end of the simulation, the sphere has almost reattained its original shape, except a spurious ridge on one side, and a dent on the other. Snapshots from the time evolution is shown in Fig. 10. A very similar end configuration which includes the spurious ridge is shown in [55], where the result has been obtained by a particle based method. The interface experiences the most stretching at the midpoint of the simulation. With both grid resolutions, the thin interface is poorly resolved. Parts of it are thus filtered out, resulting in volume loss. At 0.73% on the 50 × 50 × 50 primal grid, the volume loss is less than that obtained both with the fifth order WENO scheme and the Particle Level Set method on a 100 × 100 × 100 grid [7]. 7. Conclusive summary and future work In the present work, we have described a method for gradient augmented level set advection and reinitialization on unstructured grids. As demonstrated through the test cases in Section 6, the proposed method provides high performance on unstructured grids, both in terms of interface accuracy and volume preservation. Applied to the advection of smooth functions such as in Section 6.1, the proposed gradient limiter has little effect on the accuracy and convergence. When applied to geometrical shapes with sharp or under-resolved features, as in Sections 6.2 and 6.5, it suppresses the appearance of spurious interface artifacts. The proposed Fast Marching Method has been demonstrated to work as intended; it introduces very little numerical diffusion, and in some cases, it counteracts numerical diffusion. While the method should not be used with low-quality grids, modern meshing software can generate high-quality unstructured grids that are adequate for use with the proposed method. Note that the method is likely to perform better on a ‘native’ tetrahedral grid than the present subgrid construction, as the connectivity graph is more regular. This can be seen in Fig. 1(b), where the majority of nodes in the (2D) primal grid have six neighbors, whereas nodes with only four neighbors are evenly distributed in the subgrid. The number of neighbor nodes are important in the FMM, as discussed in Section 5.2. From the reported results, we judge the method suitable for practical flow simulations that involve a high degree of interface deformation. In extreme cases however, for instance in flows that are completely dominated by small bubbles or thin films, volume loss might become an issue. While the obvious solution is to increase the grid resolution, this is not always computationally feasible. A possible solution that has the potential to almost eliminate numerical diffusion and mass loss is the introduction of Lagrangian particles. These may track both interface location and gradient vector [9], and upon inclusion in the numerical stencils, one such particle can determine the four Bézier coefficients which are otherwise approximated by Eq. (B.4). Preliminary result have been promising, and suggest that such an approach relies much less on a dense particle distribution than e.g. the Particle Level Set method [7]. An ongoing study of the present method applied to flow simulations has shown highly encouraging results.

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

65

Acknowledgements The current work is based on a tentative implementation of a parallel semi-Lagrangian narrow-band level set method in the CFD code CDP [30]. The implementation was based on linear interpolation, both in the advection and Fast Marching Method used for reinitialization. The authors would like to express their gratitude to Associate Professor Mikael Mortensen at the University of Oslo for a critical review of the article manuscript. The current work is sponsored by The Norwegian Research Council grant NFR 1912 04/V30 ‘Wave-current-body interaction’ and the DNV Education Fund. Appendix A. Discretization by Huygens’ principle on a tetrahedron A.1. Relation to the eikonal equation Eq. (22) satisfies the eikonal equation commonly used to reinitialize the level set function:

∇Φ · ∇Φ = 1.

(A.1)

To see this, we note that ∇Φ does not change along its own direction unless colliding fronts are present. This can be seen by taking the gradient of Eq. (A.1), which yields

∇∇Φ · ∇Φ = 0.

(A.2)

As in Section 5.2, we define T to be the convex hull of the accepted simplex:

 T=

3

x∈R : x=

n +1 

ak xk ∧

k =1

n +1 

 ak = 1 ∧ ak  0 ,

(A.3)

k =1

where xk denotes a vertex coordinate, and n ∈ {0, 1, 2}. In the following, we assume n = 2, i.e., that the accepted simplex is a triangle. Assuming that Eq. (A.1) is satisfied in the whole tetrahedron, and that the gradient is continuous, the gradient through xc can be backtraced to find a point of intersection x in T that has the same gradient. For this purpose, we introduce a local Cartesian coordinate system with origin and basis vectors e 1 and e 2 in the same plane as T . Approximating the in-plane components of the gradient in T by the gradient of some interpolant P , the component along e 3 must be scaled such that the gradient has unit magnitude to satisfy Eq. (A.1). Since the close node is further from the interface than the accepted nodes, the gradient must point in the close node’s general direction, i.e., it must have a positive component along e 3 . The gradient g (x) may thus be approximated by



∂P ∂P g ( x) ≈ e1 + e2 + ∂ x1 ∂ x2

 1−

∂P ∂ x1

2

 −

∂P ∂ x2

2 e3 ,

(A.4)

assuming that the magnitude of the in-plane gradient of P is equal to or less than one. A possible numerical solution of the eikonal equation at the close node may be found by looking for a point x in T that has a gradient vector given by Eq. (A.4), whose extension from x passes through the close node xc . At such a point, an exact Taylor expansion of the distance function from this point to the close node amounts to adding the distance D = xc − x, as the gradient is of unit magnitude, and unchanging along its own direction. If f (xmin ) = P (xmin ) + D (xmin ) is a local minimum of f (see Eq. (21)) in the interior of T , and has a gradient g (xmin ) given by Eq. (A.4) which is parallel to xc − xmin , the solution to Eq. (22) complies with the eikonal equation. This is the case, and it can be seen by taking the dot product of the gradient g (xmin ) and the vector xc − xmin . At xmin , we have

xij − xmin ∂D ∂P j = = , ∂xj ∂xj D

j = 1, 2,

(A.5)

where xi is the orthogonal projection of xc onto the plane of T . The distance vector between xmin and the close node is



 min

d = xc − xmin = x1i − x1



 min

e 1 + x2i − x2



e2 + D 1 −



x1i − xmin 1 D

2

 −

x2i − xmin 2 D

2 e3 .

(A.6)

Substituting Eq. (A.5) into Eq. (A.4) and taking the dot product with Eq. (A.6) yields





g xmin · d = D .

(A.7)

Since  g  = 1 and d = D, g (xmin ) and xc − xmin must be parallel. The same analysis applies on the edges of T . Here, the gradient must be assumed to lie in the plane formed by xc and two of the accepted nodes. We previously mentioned that the gradient is unchanging along its own direction unless colliding fronts are present. If colliding fronts are present,

66

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

the minimization over T in Eq. (22) ensures that the front that emanates from xmin = xe is not blocked by other fronts emanating from T . This, of course, requires that the minimization procedure converges to the true minimum. The FMM and minimization over adjacent accepted simplices A ensure that the front emanating from xe is not blocked by fronts coming from outside of T , and thus, the gradient is continuous on the line segment between xc and xe , and the above assumptions are valid. A.2. Update procedure and discretization We will describe the update procedure assuming three accepted neighbors; the procedure is analogous for two neighbors, and for one neighbor, there is no need to search for a minimum of f , as there is only one eligible point. To solve Eq. (22) on a tetrahedron, we use the local Cartesian coordinate system with origin and basis vectors e 1 and e 2 in the plane of T . In the local coordinate system, we define the following points: xc is the close node, xi is the orthogonal projection of xc onto the plane of T , and x is some point in the plane of T , at which we hope to find a local minimum of f (x). We also define D = xc − x and h = xc − xi . P (x) is taken as the gradient limited version of Eq. (7). The setup described above is illustrated in Fig. 3. To apply Newton’s method to find the minimum of f (x), we need the first and second derivatives of P and D. For P , these are found from Eq. (5). For D, they are given by

D=

 

x1i − x1

2

2  + x2i − x2 + h2 ,

xi − x1 ∂D =− 1 , ∂ x1 D xi − x2 ∂D =− 2 , ∂ x2 D (xi − x1 )2 1 ∂2 D = − 1 3 , 2 D D ∂ x1

(x2i − x2 )2 1 ∂2 D = − , D D3 ∂ x22 (xi − x1 )(x2i − x2 ) ∂2 D =− 1 . ∂ x1 ∂ x2 D3

(A.8)

Newton’s method requires an initial guess x0 . For this, we discard the gradient information at the accepted nodes, and use a linear interpolant on the nodal function values. The gradient of this interpolant is constant everywhere, and there will thus be exactly one point where the eikonal-complying gradient, whose component along e 3 is scaled such that the gradient has unit magnitude, extends through the update node, see Section A.1. Having obtained the eikonal-complying gradient, the gradient vector is backtraced from the update node, and the point of intersection with the accepted simplex is taken as the starting point. Under certain circumstances, the linear in-plane gradient may have a magnitude greater than one, and does hence not comply with the eikonal equation. If that happens, xi is instead chosen as an initial guess. The linear in-plane gradient may be found by differentiating the linear interpolant P lin constructed by barycentric coordinates. Using Einstein’s summation convention,

P lin (λ) = λ j Φ j ,

∂ P lin 1 = T− Φ j, ji ∂ xi

(A.9)

where the matrix T is found from the triangle variant of Eq. (3)



T =

x1 1

x2 1



x3 , 1

(A.10)

x1 , x2 and x3 are the vertex positions in the 2D coordinate system formed by e 1 and e 2 , and j is a summation index over the vertices. The third gradient component is positive along e 3 , such that the linear gradient vector g lin has unit magnitude:



g

lin

−1

−1

= T j1 Φ j e 1 + T j2 Φ j e 2 +



1 1 − T− Φj j1

2

 1 2 − T− Φ j e3 , j2

(A.11)

where Einstein’s summation convention again has been invoked. Backtracing g lin from the close node to the point of intersection with the accepted triangle yields

x0 = x i −

h  g 3lin



g 1lin e 1 + g 2lin e 2 .

(A.12)

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

67

While providing an initial guess, the above equations may also be used as a standalone local, linear update rule, by updating Φ in the close node according to

      Φ xc ≈ P lin x0 + xc − x0 .

(A.13)

As mentioned above, if the square root in Eq. (A.11) is complex, x0 is replaced by xi . When this happens, however, it is an indication that the FMM has locally failed to approximate the distance function satisfactorily. When using the linear approximation above to assign values to new band members (Section 2), a complex square root is taken as a failure, and the minimum search is continued on the edges and vertices. It is interesting to note that the linear approach given by Eqs. (A.9)–(A.13) is equivalent to the first order Finite Difference approach given in [46]. In [46], a linear variation in Φ along the edges from the vertices of the accepted triangle to the close node is assumed. These provide linearly independent directions, and thus determine a constant gradient in the whole tetrahedron. Combined with the requirement of a unit magnitude gradient, the problem has two solutions; one where the gradient has a positive component along e 3 and one with a negative component, the former being the meaningful travel time estimate. These are the solutions of the quadratic equation presented in [46]. The same assumptions are made in the above equations, and empirical investigation has confirmed that the two approaches produce identical results down to machine precision. As an alternative to Eq. (A.13), it is equally acceptable to find the gradient by Eq. (A.11) and approximate Φ(xc ) by a Taylor expansion from one of the accepted nodes, which produces exactly the same result. A.3. Optimization algorithm Newton’s method for optimization is an iterative method based on a quadratic Taylor expansion around at the current iterate point, and finding the minimum of the Taylor-expansion. This yields a linear system of equations:

H i j u j = − gi ,

∂ 2 f  n x , ∂ xi ∂ x j ∂ f  n gi = x , ∂ xi

Hij =

u = xn+1 − xn ,

(A.14)

where n is the current iterate, and the summation convention is invoked. Eq. (A.14) alone may of course diverge, or converge to a local maximum or saddle point. The most common methods applied to ensure convergence to a local minimum are called line-search methods and trust region methods. TOMS611 applies the latter. Trust region methods define a ‘trusted region’, within which the objective function may adequately be approximated by a quadratic function. The region is redetermined in every iteration, based on the ratios of the objective function to the Taylor expansion. Minimization is then constrained to the trusted region. Appendix B. A Hermite interpolant for tetrahedral data To construct a tetrahedral cubic Hermite interpolant, we follow a similar procedure as with a tensor product interpolant – i.e., we build bi- and trivariate interpolants from the univariate case. B.1. The univariate Hermite interpolant We express the well-known cubic Hermite interpolant in Bernstein–Bézier form (B-form), which is convenient to work with when the number of variables increases. The univariate nth degree B-form in barycentric coordinates reads

P (λ) =



c i j B nij ,

i + j =n

B nij =

n! i! j!

j

λ1i λ2 .

(B.1)

There is only one independent variable here, as λ1 + λ2 = 1. If 0  λ1 , λ2 , the interpolant is constrained to a line segment, or a 1-simplex. The cubic Hermite interpolant, which is also the basis in the tensor product approach, is obtained by

n = 3, c 30 = f 10 , c 21 = f 10 +

1 3

(x01 − x10 ) · g 10 ,

(B.2)

68

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

where f 10 , g 10 and x10 are the function value, gradient and Cartesian position vector at λ = (λ1 , λ2 ) = (1, 0), i.e., the vertex data and coordinates of vertex 1. The remaining two coefficients – c 03 and c 12 – are found by permuting the entries of the multi-index. B.2. Extension to bivariate Hermite interpolation The bivariate B-form in barycentric coordinates reads



P (λ) =

c i jk B nijk ,

i + j +k=n

n!

B nijk =

j

i ! j !k!

λ1i λ2 λk3 ,

(B.3)

which expresses a cubic function if n = 3. If 0  λ1 , λ2 , λ3 , the interpolant is constrained to a triangle, or a 2-simplex. For λm = 0, we require that the univariate polynomial along the edge that excludes vertex m be recovered. For example, if λ3 = 0, then B nijk = B nij if k = 0, otherwise, B nijk = 0. This applies similarly if λ2 = 0, which should recover the univariate polynomial along edge 1–3. We can use this directly to determine all of the coefficients but c 111 , simply by adding a third zero-index in Eq. (6), and use index permutation to obtain the other coefficients. To find c 111 , we follow [56], where the function at the barycenter of the triangle is approximated by a Taylor expansion:

f i jk = xi jk =

1 3 1 3

( f i + f j + fk) −

1 6



 

g i · (xi − xi jk ) + g j · (x j − xi jk ) + g k · (xk − xi jk ) + O h3 ,

(xi + x j + xk ),

(B.4)

where h is a characteristic length of the triangle. The vertex data in the above equation may instead be expressed in terms of the Bézier coefficients:

f i jk ≈

1 6

(c 210 + c 201 + c 120 + c 102 + c 021 + c 012 ).

(B.5)

Equating Eq. (B.3) evaluated in λ = ( 13 , 13 , 13 ), and with n = 3 to f i jk yields

1 27



6

i + j +k=3

i ! j !k!

c i jk = f i jk .

(B.6)

Approximating f i jk by Eq. (B.5) and solving for c 111 yields

c 111 =

1 4

1

(c 210 + c 201 + c 120 + c 102 + c 021 + c 012 ) − (c 300 + c 030 + c 003 ). 6

(B.7)

It should be noted that not all of the coefficients in the tensor product approach presented in [15] are readily obtained from stored vertex data, and finite differences are used in a similar way to approximate the missing vertex data. In the present work, the requirement of compact stencils makes the order of accuracy drop by one as compared to the work of [15]. B.3. Extension to trivariate Hermite interpolation The trivariate B-form in barycentric coordinates reads

P (λ) =



c i jkl B nijkl ,

i + j +k+l=n

B nijkl =

n! i ! j !k!l!

j

λ1i λ2 λk3 λl4 ,

(B.8)

which expresses a cubic function if n = 3. If 0  λ1 , λ2 , λ3 , λ4 , the interpolant is constrained to a tetrahedron, or a 3-simplex. For λm = 0, we require that the bivariate polynomial along the face that excludes vertex m be recovered. We can use this directly to determine all of the coefficients simply by adding a fourth zero-index in Eq. (7), and use index permutation to obtain the other coefficients. The resulting interpolant is equivalent to those given in [19,56].

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

69

B.4. Gradient limiting A distance function will have discontinuous first derivatives at points that do not have a unique closest point on the interface. Representing the distance function in terms of propagating fronts (as in the FMM) these points are where two or more fronts collide, and are manifested by ‘kinks’ in the distance function. In areas of high curvature and thin films, these kinks appear near the interface. The distance function across a kink is poorly described by cubic Hermite interpolation, and the resulting overshoots tend to increase the size of small-scale features as reconstructed by the interpolant. While we can accept the degraded accuracy near kinks, we cannot accept the overshoots. Implementation of the method with a flow solver has revealed that overshoots may cause small-scale features to grow over time, which leads to results that are obviously unphysical. Instead, we prefer small scale features to be preserved, or filtered out. To accomplish this, we limit the gradients of the distance function locally. The procedure is inspired by a common formulation of cubic Hermite splines in Computer-Aided Geometric Design (CAGD). It represents a concavity/convexity constraint. First, we translate the univariate cubic Hermite interpolant to a graph in the xy plane. While this may be done by plotting y (λ1 , λ2 ) = P (λ1 , λ2 ) against x(λ1 , λ2 ) = λ1 x1 + λ2 x2 , it may also be done in a more geometric fashion by replacing the coefficients c i j by points ai j in the xy plane. This is done by setting

 ai j =

i n

x1 +

j n

 x2 , c i j ,



x(λ1 , λ2 ) =

ai j

i + j =n

n! i! j!

j

λ1i λ2 .

(B.9)

By the Binomial theorem, the x-components of the first approach are recovered. The Binomial theorem can also be used to show that

 n! j λ i λ = 1. i! j! 1 2

(B.10)

i + j =n

Since λ1 and λ2 are always positive between x1 and x2 , any point on the graph between x1 and x2 lies within the convex hull of the points ai j . The x-components of ai j are evenly spaced along the line segment between x1 and x2 , with spacing |x2 −x1 | 3

. The y-components c i j have a very useful geometric interpretation. For the cubic Hermite polynomial, we have

c 30 = f (x1 ), c 03 = f (x2 ), c 21 − c 30 df 1 (x 3 2

− x1 )

c 03 − c 12 1 (x − x1 ) 3 2

= =

dx df dx

(x1 ), (x2 ).

(B.11)

The piecewise linear graph through a30 , a21 , a12 and a03 is the control polygon of the cubic Hermite interpolant expressed as a Bézier curve – a polygon describing the general shape of the curve. The points are referred to as control points. From Eq. (B.11), it can be seen that this polygon has end-point function values and first derivatives coinciding with the curve. A graphical explanation of why the cubic Hermite interpolant often over-predicts the size of small-scale features can be seen in Fig. 11(a). Here, the function and gradient of a 1D distance function are sampled near a small-scale feature – say, a 1D bubble. Due to the small size, the distance function exhibits a kink near the interface. To reproduce the sampled gradient, a21 must lie outside the actual distance function’s graph. The convex hull of the points ai j thus encompasses an area outside the distance function’s graph – as does the interpolant. As can be seen, this happens if the distance function and gradient are sampled closer than 13 h from a kink, where h is the grid spacing. In Fig. 11(a), this gives the interpolant a zero-contour outside the original distance function, and the size of the feature is thus overestimated. One possible way to remedy the overshoot is to decrease the y-coordinate of a21 so that it lies on the extension of the line segment between a12 and a03 : c 21 = 2c 12 − c 03 . Then, a21 lies on the distance function’s graph, see Fig. 11(b). This represents one end of the bounds for c 21 to preserve convexity/concavity; the other is c 21 = 12 (c 30 + c 12 ), which we do not use. To identify problematic end-point gradients that require limiting, we look for intersections between the extension of one end-segment of the control polygon with the other. Such an intersection is identified by looking at differences in slopes between the control points, represented by D 1 to D 5 in Eq. (9). We will now demonstrate the mechanisms of the gradient limiter by a simple experiment, which is illustrated in Fig. 12. Starting with a convex control polygon, Fig. 12(a), and gradually increasing the left end-point gradient, the gradient limiter will start limiting the right end-point gradient when the control polygon is no longer convex to retain convexity, Fig. 12(b). Continuing increasing the left end-point gradient, both end-point gradients of the gradient limited control polygon will eventually equal the linear gradient between the end-points. Further increasing the left end-point gradient will result in an S-shaped control polygon. At this point, an edge

70

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

Fig. 11. Elimination of overshoot by gradient limiting.

Fig. 12. Mechanisms of the gradient limiter.

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

71

may potentially have three interface intersections, which is a level of subgrid detail beyond what may physically be resolved in a satisfactory manner. We thus limit both end-point gradients to the linear gradient, Fig. 12(c) which is taken care of in the last two statements of Algorithm 1. Algorithm 1 is designed for the univariate cubic Hermite interpolant. But as the multivariate versions reduce to the univariate interpolant along edges, Algorithm 1 may be used edge-wise on the multivariate interpolants. This will ensure convexity/concavity along all edges. References [1] C.W. Hirt, B.D. Nichols, Volume of fluid (vof) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1) (1981) 201–225. [2] M. Huang, B. Chen, L. Wu, A slic–vof method based on unstructured grid, Microgravity Sci. Technol. 22 (3) (2010) 305–314. [3] X. Zheng, J. Lowengrub, A. Anderson, V. Cristini, Adaptive unstructured volume remeshing – ii: Application to two- and three-dimensional level-set simulations of multiphase flow, J. Comput. Phys. 208 (2) (2005) 626–650. [4] M. Herrmann, Refined level set grid method for tracking interfaces, in: Annual Research Briefs, Center for Turbulence Research, NASA Ames/Stanford University, 2005, pp. 3–18. [5] X. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1) (1994) 200–212. [6] M. Dumbser, M. Käser, Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems, J. Comput. Phys. 221 (2) (2007) 693–723. [7] D. Enright, R. Fedkiw, J. Ferziger, I. Mitchell, A hybrid particle level set method for improved interface capturing, J. Comput. Phys. 183 (1) (2002) 83–116. [8] D. Enright, F. Losasso, R. Fedkiw, A fast and accurate semi-Lagrangian particle level set method, Comput. Struct. 83 (6) (2005) 479–490. [9] S. Ianniello, A. Di Mascio, A self-adaptive oriented particles Level-Set method for tracking interfaces, J. Comput. Phys. 229 (4) (2010) 1353–1380. [10] E. Olsson, G. Kreiss, A conservative level set method for two phase flow, J. Comput. Phys. 210 (1) (2005) 225–246. [11] O. Desjardins, V. Moureau, H. Pitsch, An accurate conservative level set/ghost fluid method for simulating turbulent atomization, J. Comput. Phys. 227 (18) (2008) 8395–8416. [12] Y. Chang, T. Hou, B. Merriman, S. Osher, A level set formulation of Eulerian interface capturing methods for incompressible fluid flows, J. Comput. Phys. 124 (2) (1996) 449–464. [13] H. Takewaki, A. Nishiguchi, T. Yabe, Cubic interpolated pseudo-particle method (cip) for solving hyperbolic-type equations, J. Comput. Phys. 61 (2) (1985) 261–268. [14] T. Yabe, F. Xiao, T. Utsumi, The constrained interpolation profile method for multiphase analysis, J. Comput. Phys. 169 (2) (2001) 556–593. [15] J. Nave, R. Rosales, B. Seibold, A gradient-augmented level set method with an optimally local, coherent advection scheme, J. Comput. Phys. 229 (10) (2010) 3802–3827. [16] M. Raessi, J. Mostaghimi, M. Bussmann, Advecting normal vectors: A new method for calculating interface normals and curvatures when modeling two-phase flows, J. Comput. Phys. 226 (1) (2007) 774–797. [17] D. Salac, The augmented fast marching method for level set reinitialization, Bulletin of the American Physical Society 56. [18] N. Tanaka, T. Yamasaki, T. Taguchi, Accurate and robust fluid analysis using cubic interpolation with volume/area coordinates (civa) method on unstructured grids, JSME Int. J. Ser. B 47 (4) (2004) 672–680. [19] N. Tanaka, Development of a highly accurate interpolation method for mesh-free flow simulations i. Integration of gridless, particle and cip methods, Int. J. Numer. Methods Fluids 30 (8) (1999) 957–976. [20] O. Desjardins, H. Pitsch, A spectrally refined interface approach for simulating multiphase flows, J. Comput. Phys. 228 (5) (2009) 1658–1677. [21] S. Ii, M. Shimuta, F. Xiao, A 4th-order and single-cell-based advection scheme on unstructured grids using multi-moments, Comput. Phys. Commun. 173 (1–2) (2005) 17–33. [22] M. Lentine, J.T. Grétarsson, R. Fedkiw, An unconditionally stable fully conservative semi-Lagrangian method, J. Comput. Phys. 230 (8) (2011) 2857–2879. [23] J. Tsitsiklis, Efficient algorithms for globally optimal trajectories, IEEE Trans. Autom. Control 40 (9) (1995) 1528–1538. [24] J. Sethian, A. Vladimirsky, Ordered upwind methods for static Hamilton–Jacobi equations: Theory and algorithms, SIAM J. Numer. Anal. 41 (1) (2003) 325–363. [25] J. Grooss, J. Hesthaven, A level set discontinuous Galerkin method for free surface flows, Comput. Methods Appl. Mech. Eng. 195 (25) (2006) 3406–3429. [26] E. Marchandise, J. Remacle, N. Chevaugeon, A quadrature-free discontinuous Galerkin method for the level set equation, J. Comput. Phys. 212 (1) (2006) 338–357. [27] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49. [28] D. Adalsteinsson, J. Sethian, A fast level set method for propagating interfaces, J. Comput. Phys. 118 (2) (1995) 269–277. [29] P. Gomez, J. Hernandez, J. Lopez, On the reinitialization procedure in a narrow-band locally refined level set method for interfacial flows, Int. J. Numer. Methods Eng. 63 (10) (2005) 1478–1512. [30] F. Ham, K. Mattsson, G. Iaccarino, Accurate and stable finite volume operators for unstructured flow solvers, in: Annual Research Briefs, Center for Turbulence Research, NASA Ames/Stanford University, 2006, pp. 243–261. [31] F. Ham, K. Mattsson, G. Iaccarino, P. Moin, Towards time-stable and accurate les on unstructured grids, in: Complex Effects in Large Eddy Simulations, 2007, pp. 235–249. [32] M. Lai, L. Schumaker, Spline Functions on Triangulations, vol. 110, Cambridge University Press, 2007. [33] C. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys. 77 (2) (1988) 439–471. [34] J. Sethian, Fast marching methods, SIAM Rev. (1999) 199–235. [35] M. Crandall, P. Lions, Viscosity solutions of Hamilton–Jacobi equations, Trans. Am. Math. Soc. (1983) 1–42. [36] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1) (1994) 146–159. [37] E. Dijkstra, A note on two problems in connexion with graphs, Numer. Math. 1 (1) (1959) 269–271. [38] R. Segdewick, K. Wayne, Algorithms, 4th edition, Addison-Wesley Professional, 2011. [39] C. Huygens, Traité de la lumière [1690], Project Gutenberg, 2005. [40] H. Piaggio, The mathematical theory of Huygens’ principle, Nature 145 (1940) 531–532. [41] S. Fomel, A variational formulation of the fast marching eikonal solver, in: SEP-95: Stanford Exploration Project, 1997, pp. 127–147. [42] P.G. Lelièvre, C.G. Farquharson, C.A. Hurich, Computing first-arrival seismic traveltimes on unstructured 3-d tetrahedral grids using the fast marching method, Geophys. J. Int. 184 (2) (2011) 885–896. [43] A. de Witte, Equivalence of Huygens’ principle and Fermat’s principle in ray geometry, Am. J. Phys. 27 (1959) 293–301. [44] J. Sethian, Evolution, implementation, and application of level set and fast marching methods for advancing fronts, J. Comput. Phys. 169 (2) (2001) 503–555.

72

A. Bøckmann, M. Vartdal / Journal of Computational Physics 258 (2014) 47–72

[45] D.M. Gay, Algorithm 611: Subroutines for unconstrained minimization using a model/trust-region approach, ACM Trans. Math. Softw. 9 (4) (1983) 503–524. [46] J. Sethian, A. Vladimirsky, Fast methods for the eikonal and related Hamilton–Jacobi equations on unstructured meshes, Proc. Natl. Acad. Sci. 97 (11) (2000) 5699. [47] D. Chopp, Some improvements of the fast marching method, SIAM J. Sci. Comput. 23 (1) (2001) 230–244. [48] M. Herrmann, A domain decomposition parallelization of the fast marching method, in: Annual Research Briefs, Center for Turbulence Research, NASA Ames/Stanford University, 2003, pp. 213–225. [49] M. Sussman, E. Fatemi, An efficient, interface-preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow, SIAM J. Sci. Comput. 20 (4) (1999) 1165–1191. [50] J. Bell, P. Colella, H. Glaz, A second-order projection method for the incompressible Navier–Stokes equations, J. Comput. Phys. 85 (2) (1989) 257–283. [51] R. LeVeque, High-resolution conservative algorithms for advection in incompressible flow, SIAM J. Numer. Anal. (1996) 627–665. [52] S. Zalesak, Fully multidimensional flux-corrected transport algorithms for fluids, J. Comput. Phys. 31 (3) (1979) 335–362. [53] M. Herrmann, A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids, J. Comput. Phys. 227 (4) (2008) 2674–2706. [54] P.A. Gremaud, C.M. Kuster, Z. Li, A study of numerical methods for the level set approach, Appl. Numer. Math. 57 (5–7) (2007) 837–846. [55] S. Hieber, P. Petros Koumoutsakos, A Lagrangian particle level set method, J. Comput. Phys. 210 (1) (2005) 342–367. [56] P. Ciarlet, The Finite Element Method for Elliptic Problems, vol. 4, North Holland, 1978.