Extending fields in a level set method by solving a biharmonic equation

TL;DR

This paper introduces a biharmonic equation-based method for extending fields in level set methods, which implicitly handles interfaces and ensures smooth, accurate extensions with simple implementation and symmetry preservation.

Contribution

It presents a novel biharmonic extension approach that implicitly manages interfaces and improves smoothness and accuracy in level set methods.

Findings

Effective at producing smooth, accurate extensions near interfaces

Handles symmetry and periodicity naturally

Simple implementation with efficient solvers

Abstract

We present an approach for computing extensions of velocities or other fields in level set methods by solving a biharmonic equation. The approach differs from other commonly used approaches to velocity extension because it deals with the interface fully implicitly through the level set function. No explicit properties of the interface, such as its location or the velocity on the interface, are required in computing the extension. These features lead to a particularly simple implementation using either a sparse direct solver or a matrix-free conjugate gradient solver. Furthermore, we propose a fast Poisson preconditioner that can be used to accelerate the convergence of the latter. We demonstrate the biharmonic extension on a number of test problems that serve to illustrate its effectiveness at producing smooth and accurate extensions near interfaces. A further feature of the method is the natural way in which it deals with symmetry and periodicity, ensuring through its construction that the extension field also respects these symmetries.

Tables

Table 1. Table 1: Iteration counts, errors and estimated orders of convergence for the biharmonic extension over a sequence of grids for Example 1. Dirichlet, Neumann and mixed boundary conditions are tested. All errors are calculated with respect to the reference field ( 13 ) in a neighbourhood of 4 grid cells of the interface.

Grid	$128\times 128$	$256\times 256$	$512\times 512$	$1024\times 1024$
	Biharmonic extension (Dirichlet)
Iterations	50	112	238	543
Error	6.28e-02	1.77e-02	4.71e-03	1.18e-03
Est. order	–	1.83	1.91	2.00
	Biharmonic extension (Neumann)
Iterations	50	117	253	558
Error	5.70e-02	1.57e-02	4.10e-03	1.02e-03
Est. order	–	1.86	1.94	2.00
	Biharmonic extension (mixed)
Iterations	53	120	257	579
Error	5.57e-02	1.55e-02	4.09e-03	1.02e-03
Est. order	–	1.84	1.92	2.00

Table 2. Table 2: Runtimes, condition numbers, errors and estimated orders of convergence over a sequence of grids for Example 1 for the biharmonic extension (Neumann conditions) and the constant, linear and quadratic extrapolations. All errors are calculated with respect to the reference field ( 13 ) in a neighbourhood of 4 grid cells of the interface.

Grid	$128\times 128$	$256\times 256$	$512\times 512$	$1024\times 1024$
	Biharmonic extension (Neumann)
Runtime (direct)	0.1	0.4	2.3	12.4
Runtime (iterative)	0.1	0.9	13.2	131.0
Cond. Est.	2.4e+07	3.6e+08	5.7e+09	9.1e+11
Error	5.70e-02	1.57e-02	4.10e-03	1.02e-03
Est. Order	–	1.86	1.94	2.00
	Constant extrapolation
Runtime	0.2	1.0	13.0	128.5
Error	1.10e-01	6.06e-02	3.14e-02	1.57e-02
Est. Order	–	0.86	0.95	1.00
	Linear extrapolation
Runtime	0.3	2.0	25.8	258.9
Error	2.09e-02	5.21e-03	1.34e-03	3.38e-04
Est. Order	–	2.00	1.96	1.99
	Quadratic extrapolation
Runtime	0.5	3.0	40.8	392.5
Error	3.19e-03	4.58e-04	5.41e-05	6.73e-06
Est. Order	–	2.80	3.08	3.01

Table 3. Table 3: Runtimes, errors and estimated orders of convergence over a sequence of grids for Example 2 for the biharmonic extension (periodic conditions) and the constant, linear and quadratic extrapolations. All errors are calculated with respect to the reference field ( 14 ) in a neighbourhood of 4 grid cells of the interface.

Grid	${32}^3$	${64}^3$	${128}^3$	${256}^3$
	Biharmonic extension (Periodic)
Iterations	21	46	112	276
Runtime	0.1	1.5	37.6	733.9
Error	1.08e-01	3.42e-02	1.07e-02	2.84e-03
Est. Order	–	1.66	1.68	1.91
	Constant extrapolation
Runtime	0.2	5.3	85.9	1269.9
Error	1.29e-01	6.87e-02	3.55e-02	1.77e-02
Est. Order	–	0.91	0.95	1.00
	Linear extrapolation
Runtime	0.4	10.3	150.4	2258.8
Error	5.17e-02	1.14e-02	2.93e-03	7.42e-04
Est. Order	–	2.18	1.96	1.98
	Quadratic extrapolation
Runtime	0.5	15.6	229.0	3420.3
Error	2.99e-02	3.31e-03	4.52e-04	7.66e-05
Est. Order	–	3.18	2.87	2.56

Table 4. Table 4: Iteration counts, errors and estimated orders of convergence over a sequence of grids for Example 3. Dirichlet and Neumann boundary conditions are tested, and errors are reported for inside and outside the annulus. All errors are calculated with respect to the reference field ( 15 ) in a neighbourhood of 4 grid cells of the interface.

Grid	$128\times 128$	$256\times 256$	$512\times 512$	$1024\times 1024$
	Dirichlet
Iterations	56	123	260	529
Error outside	6.15e-02	1.74e-02	4.66e-03	1.22e-03
Est. Order outside	–	1.83	1.90	1.94
Error inside	9.73e-02	2.80e-02	7.76e-03	2.07e-03
Est. Order inside	–	1.80	1.85	1.90
	Neumann
Iterations	57	124	256	522
Error outside	4.06e-03	1.07e-03	3.00e-04	7.91e-05
Est. Order outside	–	1.92	1.83	1.92
Error inside	9.73e-02	2.80e-02	7.76e-03	2.07e-03
Est. Order inside	–	1.80	1.85	1.90

Table 5. Table 5: Runtimes, errors and ratios of errors over a sequence of grids for Example 2. Dirichlet and Neumann boundary conditions are tested, and errors and ratios are reported for inside and outside the annulus. All errors are calculated with respect to the reference field ( 15 ) in a neighbourhood of 4 grid cells of the interface.

Grid	$128\times 128$	$256\times 256$	$512\times 512$	$1024\times 1024$
	Biharmonic extension (Neumann)
Runtime (direct)	0.1	0.4	2.1	11.3
Runtime (iterative)	0.2	0.9	13.1	121.4
Cond. Est.	2.0e+07	3.1e+08	4.9e+09	7.8e+10
Error outside	4.06e-03	1.07e-03	3.00e-04	7.91e-05
Est. Order outside	–	1.92	1.83	1.92
Error inside	9.73e-02	2.80e-02	7.76e-03	2.07e-03
Est. Order inside	–	1.80	1.85	1.90
	Constant extrapolation
Runtime	0.2	1.0	12.9	128.3
Error outside	4.48e-02	2.33e-02	1.18e-02	6.05e-03
Est. Order outside	–	0.94	0.99	0.96
Error inside	4.92e-02	2.45e-02	1.26e-02	6.75e-03
Est. Order inside	–	1.00	0.96	0.90
	Linear extrapolation
Runtime	0.3	2.0	25.7	259.7
Error outside	7.64e-04	1.77e-04	4.26e-05	1.06e-05
Est. Order outside	–	2.11	2.06	2.00
Error inside	1.06e-03	2.27e-04	5.65e-05	1.62e-05
Est. Order inside	–	2.22	2.01	1.81
	Quadratic extrapolation
Runtime	0.5	3.0	39.4	389.5
Error outside	1.43e-04	1.83e-05	2.40e-06	3.05e-07
Est. Order outside	–	2.96	2.93	2.97
Error inside	1.49e-03	1.90e-04	2.21e-05	2.77e-06
Est. Order inside	–	2.97	3.10	3.00

Table 6. Table 6: Iteration counts, condition numbers, errors and estimated orders of convergence over a sequence of grids for Example 4 with periodic conditions on the top and bottom. Two interfaces are tested: the irregular interface pictured in Figure 7 consisting of 10 terms of scaled and shifted sine waves with increasing wavenumber, and a simpler interface consisting of just the first term. All errors are calculated with respect to the reference field ( 16 ) in a neighbourhood of 4 grid cells of the interface.

Grid	$192\times 64$	$384\times 128$	$768\times 256$	$1536\times 512$
	10 terms
Iterations	86	185	424	979
Runtime (iterative)	0.3	1.1	19.6	200.3
Runtime (direct)	0.0	0.2	1.4	7.0
Cond. Est.	5.8e+08	9.1e+09	1.4e+11	2.3e+12
Error	2.45e-02	9.27e-03	2.96e-03	1.03e-03
Est. Order	–	1.40	1.65	1.52
	1 term
Iterations	38	70	131	269
Runtime (iterative)	0.5	0.5	5.9	52.3
Error	1.54e-02	4.66e-03	1.30e-03	3.52e-04
Error Order	–	1.72	1.84	1.88

Table 7. Table 7: Iteration counts, condition numbers, errors and estimated orders of convergence over a sequence of grids for Example 4 with Neumann conditions on the top and bottom. Two interfaces are tested: the irregular interface pictured in Figure 8 consisting of 10 terms of scaled and shifted sine waves with increasing wavenumber, and a simpler interface consisting of just the first term. All errors are calculated with respect to the reference field ( 17 ) in a neighbourhood of 4 grid cells of the interface.

Grid	$192\times 64$	$384\times 128$	$768\times 256$	$1536\times 512$
	10 terms
Iterations	97	212	496	1176
Runtime (iterative)	0.3	1.3	17.9	206.6
Runtime (direct)	0.1	0.2	1.1	6.3
Cond. Est.	5.6e+08	9.1e+09	1.5e+11	2.3e+12
Error	8.82e-03	3.11e-03	1.11e-03	3.76e-04
Est. Order	–	1.50	1.48	1.57
	1 term
Iterations	43	88	180	384
Runtime (iterative)	0.1	0.7	7.2	68.0
Error	5.88e-03	1.63e-03	4.54e-04	1.32e-04
Error Order	–	1.85	1.84	1.79

Equations

$$V_n=\mathbf{F}\cdot \mathbf{n}\text{on }\partial \Omega .$$

$$\frac{\partial \varphi }{\partial t}+\mathbf{F}\cdot \nabla \varphi =0.$$

$$\nabla f\cdot \nabla \varphi =0$$

$$\frac{\partial f}{\partial \tau }+H(\varphi )\nabla f\cdot \mathbf{n}=0$$

$$f({\mathbf{x}}_k)=f_k,\text{for }k=1,\dots ,N.$$

$$f(\mathbf{x})=\sum _{k=1}^N{c}_k\Phi (\Vert \mathbf{x}-{\mathbf{x}}_k\Vert ),$$

$$\Phi (r)=r^2\log (r)$$

$${\iint }_{{\mathbb{R}}^2}(f_{xx}^2+2f_{xy}^2+f_{yy}^2)\text{d}x\text{d}y$$

$${\nabla }^{4}f=0\text{in }D\backslash \{{\mathbf{x}}_1,\dots ,{\mathbf{x}}_N\}$$

$$f=0\text{and}{\nabla }^{2}f=0\text{on }\partial D$$

$$f({\mathbf{x}}_k)=f_k,\text{for }k=1,\dots ,N.$$

$$f=f_k\text{and}\partial f/\partial n=0\text{on }{\Gamma }_k$$

$$20f_{i,j}-8f_{i-1,j}-8f_{i,j+1}-8f_{i,j-1}-8f_{i+1,j}+2f_{i+1,j+1}+2f_{i-1,j-1}$$

$$+2f_{i+1,j-1}+2f_{i-1,j+1}+f_{i+2,j}+f_{i-2,j}+f_{i,j+2}+f_{i,j-2}=0$$

$$f(x,y)=\cos (x)\sin (y)$$

$$\Omega =\left\{(x,y)\in D|\min \left(\sqrt{(x-0.8{)}^2+y^2}-1,\sqrt{(x+0.8{)}^2+y^2}-1\right)<0\right\}.$$

$$f(x,y,z)=\cos (x)\sin (y)\sin (\pi /4-z)$$

$$\Omega =\left\{(x,y,z)\in D|\min \left(\sqrt{(x-0.8{)}^2+y^2+z^2}-1,\sqrt{(x+0.8{)}^2+y^2+z^2}-1\right)<0\right\}.$$

$$f(x,y)=y/\log (1+\sqrt{x^2+y^2})$$

$$\Omega =\left\{(x,y)\in D|1/2<\sqrt{x^2+y^2}<1\right\}.$$

$$f(x,y)=x+\cos (2\pi (y-0.25))/5$$

$$f(x,y)=x+\cos (\pi y)/5$$

$$\frac{\partial T}{\partial t}={\nabla }^{2}T\text{in}D.$$

$$V_n=-\left(\frac{\partial T_{\ell }}{\partial n}-\frac{\partial T_s}{\partial n}\right)=-\left[\frac{\partial T}{\partial n}\right]\text{on}\partial \Omega .$$

$$T=-\sigma \kappa \text{on}\partial \Omega ,$$

$$T\left(\mathbf{x},0\right)=\begin{dcases*}0,&$\mathbf{x}\in\Omega$\\ -\frac{1}{\beta},&$\mathbf{x}\in D\backslash\Omega$\end{dcases*},$$

Full text

We demonstrate the biharmonic extension on a number of test problems that serve to illustrate its effectiveness at producing smooth and accurate extensions near interfaces. A further feature of the method is the natural way in which it deals with symmetry and periodicity, ensuring through its construction that the extension field also respects these symmetries.

††journal: Journal of Computational Physics

1 Introduction

Moving boundary problems, or interface problems, arise in the mathematical modelling of many interesting natural phenomena, including Hele-Shaw flows [1, 2, 3, 4], multiphase fluid flows [5], thin film flows [6, 7, 8], melting and solidification problems [9, 10, 11, 12], melanoma growth [13] and many others.

Schematics of two particular configurations of a moving boundary problem on a two-dimensional domain $D$ are shown in Figure 1. We suppose that $D$ is partitioned into two regions, and that the interface $\partial\Omega$ between the two evolves under a normal velocity given by

[TABLE]

In the shaded region $\Omega$ the velocity field $\mathbf{F}$ is calculated from the solution of a field equation, which we suppose is a parabolic or elliptic partial differential equation. (We do not consider hyperbolic equations in this paper, since as noted by Gibou et al. [14] they require further considerations to ensure that conservation properties are conserved, in order that the correct speed of propagation is calculated in the presence of shocks.) Having calculated the velocity field $\mathbf{F}$, the normal velocity of the interface $V_{n}$ is then determined from (1). In this way, the interface evolves over time and, depending on the phenomena under consideration, may develop fingering instabilities, splitting, merging, and other interesting features. Accurately tracking the evolution of the interface over time is a key requirement of any numerical solution method for such a problem.

The level set method is one numerical method that has proved very popular for solving interface problems [15, 16, 17]. Rather than developing an explicit representation for the interface location, the level set method represents the interface implicitly, as the zero level set of an auxiliary function $\phi$ – the so-called level set function. This approach has the distinct advantage that the location of the interface, as well as splitting, merging, and other such complex behaviours as the interface evolves, are all handled naturally by the level set function discretised on a standard rectangular grid. The evolution of $\phi$ under the external velocity field $\mathbf{F}$ is governed by the level set equation

[TABLE]

A typical numerical scheme proceeds by solving the field equation to determine the velocity field, followed by updating the level set according to (2), discretised using finite differences. There is, however, a catch. Although the interfacial velocity $V_{n}$ has physical significance only on the interface itself, the level set equation (2) requires a field $\mathbf{F}$ defined on all $D$. Hence, a means of extrapolating, or extending, the velocity field $\mathbf{F}$ from $\Omega$ to $D$ is required to complete the level set implementation. This extension velocity need not have any physical significance itself, and is required only to be a smooth, well-defined field on all of $D$ that accurately represents the true velocity in a small neighbourhood of the interface.

Many numerical schemes have been proposed for constructing extension velocities. One approach is to solve the equation

[TABLE]

by a process known as the Fast Marching Method [18], which sweeps forward in the upwind direction from the interface, building the extension within an evolving narrow band. This approach is motivated by the desire to maintain the level set $\phi$ as a signed distance function, satisfying $\|\nabla\phi\|=1$, under evolution by (2). Higher order variants of the method can also be constructed provided appropriate care is taken with the direction of propagation near the interface for the higher order normal derivatives [19].

An alternative approach to constructing extension velocities is to solve an advection equation

[TABLE]

to steady-state in pseudo-time $\tau$, where $H$ denotes the Heaviside function [20]. The resulting field is constant along the normal directions, carrying the interfacial value along these characteristics. Higher order schemes can be built upon this foundation by a sequence of extrapolations, first with higher order normal derivatives, through lower order normal derivatives and finally the value itself [21]. Alternatively, the steady version of (3) may be solved directly using a fast sweeping method that carefully chooses and adjusts the required upwind stencils and systematically alternates through different node orderings to cover all characteristic directions [22].

A common feature of these approaches is the need to determine certain explicit properties related to the interface, such as its location, the location of the nearest node to it, or the value of the field on the interface.

This is in contrast with the level set method itself, which deals with the interface in an entirely implicit fashion. Hence we are motivated to pursue a velocity extension approach that also deals with the interface implicitly, and does not rely on an explicit representation of the interface.

A further point we address concerns problems possessing symmetry or periodicity. For these problems, it is natural to exploit this symmetry in the implementation of the level set method by solving over a reduced computational domain. This approach is valid provided the extension field $\mathbf{F}$ also possesses the required symmetry. Hence, it is an important consideration when constructing velocity extensions for such problems that they conform to the prescribed symmetry or periodicity conditions on the domain boundary.

In this paper we present an extension method that deals only with the interface implicitly, making for a scheme that is straightforward to implement and fully in keeping with the spirit of the level set method. Furthermore, the extension field naturally conforms to any desired symmetry or periodicity conditions on the boundary of the computational domain. The method should be useful for extending velocity fields, but also temperatures, mass fractions, or other fields relevant to the application. In the next section we outline the method, introducing the theoretical background, and we also propose a simple numerical implementation. In Section 3 we present five examples that illustrate the effectiveness of the method. Finally, we draw our conclusions in Section 4.

2 Biharmonic extension

2.1 Formulation of our approach

In the problem as described thus far, values of some vector field $\mathbf{F}$ are known on one side of an interface and must be extended in a sensible way to the entire domain. From this point onwards, we consider only a scalar field, which may be a component of $\mathbf{F}$, or more generally any scalar field for which an extension is required. Values of this field are known at all nodes of the finite difference grid in $\Omega$. The extension problem is then: Given known values $f_{k}$ for $k=1,\ldots,N$ at nodes $\mathbf{x}_{k}\in\Omega$, construct a suitably smooth function $f:D\to\mathbb{R}$ such that

[TABLE]

To motivate our approach, we begin by considering the simpler problem of interpolation, such as would be the case if the domain was partitioned as in Figure 1(b), that is, with the field values $f_{k}$ known on the outside of $\partial\Omega$. The problem of determining field values on the inside is then an interpolation, since the interpolation points lie within the convex hull of the prescribed data. In two dimensions, a natural choice of interpolant is the thin plate spline [23]. Formulated in terms of radial basis functions (RBFs) [24], the thin plate spline interpolant takes the form

[TABLE]

where $\Phi:\mathbb{R}^{+}\to\mathbb{R}$, a function of the distance between points, is given by

[TABLE]

and the $N$ coefficients $c_{k}$ are determined by enforcing the $N$ interpolation constraints (4).

The terminology “thin plate spline” derives from one application of this interpolant, which when augmented by an affine polynomial term (i.e. $f(\mathbf{x})+a_{0}+b_{1}x+b_{2}y$) describes the shape of a thin metal plate, infinite in extent, which is constrained to pass through the set of interpolation points. In particular, the function minimises the “bending energy”

[TABLE]

over all possible interpolants $f$ [23].

The thin plate spline interpolant can in principle be determined directly by solving an $N\times N$ linear system for the coefficients $c_{k}$ in (5) (or, if the affine polynomial term is included, an $(N+3)\times(N+3)$ linear system, which allows for the imposition of three further constraints that force the solution to be asymptotically flat as $\|\mathbf{x}\|\to\infty$). There are, however, several limitations that render this approach impractical for the present application of computing extension velocities. First, it is well known that the conditioning of the RBF linear system is typically very poor, and increasingly so for larger $N$. Second, whenever the RBF has global support (as it does for the thin plate spline), the coefficient matrix is dense, and hence storage and factorisation costs are prohibitive for large $N$. Third, in general, computing an extension velocity requires extrapolation rather than interpolation, such as when the domain is partitioned as in Figure 1(a), and the thin plate spline is unbounded away from the prescribed data.

These limitations may be overcome by considering the thin plate spline interpolant from an alternative point of view [25]. The thin plate spline RBF (6) is (to within a multiplicative constant) a fundamental solution of the biharmonic equation $\nabla^{4}\Phi=0$ in $\mathbb{R}^{2}$. Hence, by superposition, the linear combination (5) also solves the biharmonic equation in $\mathbb{R}^{2}$, except possibly at the points $\mathbf{x}_{k}$ themselves. Being a function defined on all of $\mathbb{R}^{2}$, it satisfies no particular boundary conditions and indeed as already noted, it is not bounded for large $\|\mathbf{x}\|$. However, if we instead consider the problem on the finite domain $D\subset\mathbb{R}^{2}$, we may now enforce boundary conditions on $f$, which can be used to “clamp” the thin plate spline, for example by enforcing that it vanishes on $\partial D$.

The idea then, is to determine the extension field $f$ as the solution to the biharmonic equation on $D\backslash\{\mathbf{x}_{1},\ldots,\mathbf{x}_{N}\}$, subject to suitable boundary conditions on $\partial D$ and also subject to the interpolation constraints (4). We shall call this function $f$ the biharmonic extension which, supposing for the present time, that homogeneous Dirichlet boundary conditions are imposed, satisfies

[TABLE]

subject to

[TABLE]

and

[TABLE]

We note that this fourth order problem is well-posed with the specification of both (9) and (10). Indeed, the problem could instead be posed by surrounding each interpolation point $\mathbf{x}_{k}$ with a circle $\Gamma_{k}$ of radius $\epsilon$, and imposing

[TABLE]

instead of (10). The unique biharmonic extension would then be obtained in the limit $\epsilon\to 0$ [25].

Solving (8) – (10) to determine the biharmonic extension is very convenient in the context of the level set method. The solution can be computed on the same finite difference grid already in place for the level set equation (2). Just as with the level set equation, there is no need to take special care with finite difference stencils near the interface. Stencils can, and do, cross over the interface, and indeed this is the whole point of computing the extension – it carries information from one side of the interface to the other.

Numerically, we simply solve the discrete biharmonic equation

[TABLE]

at every node $(x_{i},y_{j})$ in $D\backslash\Omega$. Stencils near the domain boundary are modified according to the boundary conditions (9). Where the stencil crosses the interface, the referenced values on the other side are known: they are the values $f_{k}$ that we wish to extend. In this way, the biharmonic extension naturally incorporates two levels of information on one side of the interface (the values $f_{i-1,j}$ and $f_{i-2,j}$, say) and constructs a smooth extension that follows the trend of the data across the interface to the other side. The resulting extension is, by construction, continuous in both value and derivative.

We also emphasise that at no time is the location of the interface, nor the value of the field on the interface, required in constructing the extension. The biharmonic extension deals only with the field values on the level set grid, eliminating any need to explicitly locate, or map values to, the interface before extending across to the other side. Indeed, a significant advantage of the biharmonic extension is its simplicity of implementation, as we now discuss.

2.2 Efficient computation of biharmonic extension

As already mentioned, forming the extension field by solving the biharmonic equation rather than with the RBF equation (5) avoids the dense linear systems associated with RBF-based interpolation. Instead, the coefficient matrix is the discrete biharmonic matrix arising from imposing (2.1) at each node. Nevertheless, it is important to consider how this sparse, symmetric positive definite linear system can be solved in order to build the extension field efficiently.

For problems where the number of nodes at which the extension is required is small compared to the total problem size (i.e. where $\Omega$ comprises most of the computational domain $D$), it is efficient to simply form and factorise the coefficient matrix using standard sparse algorithms (e.g. MATLAB’s sparse chol command). For extensions involving a large number of nodes, the matrix fill-in incurred by factorisation may render this approach unviable, and alternative solution methods may be required. Fast finite difference methods for biharmonic equations over irregular domains have received some attention in the literature, for example [26, 27, 28]. The idea is to embed the irregular domain in a rectangular region and apply fast Poisson solvers over the rectangle (exploiting the fact that the biharmonic operator is the square of the Laplacian), with methods differing in the way the embedding is carried out. Other approaches for building thin plate spline interpolants through the biharmonic equation have used adaptive quadtree and octtree grids [25]. Further ideas can be found from the field of image “inpainting”, where Poisson or biharmonic equations are often solved to “fill in” missing details from images. Here, methods include multigrid methods [29] and methods based on conjugate gradients and pseudo-timestepping [30]. Hence, there is great scope for further research into developing fast methods for computing these extension fields.

In the present context we have an appropriate rectangular region $D$, and so we opt for a method based on a simple embedding coupled with a fast Poisson solver. The mapping from $D\backslash\Omega$ to $D$ is achieved by augmenting the right hand side vector of the linear system with zeros, while the mapping back to $D\backslash\Omega$ simply omits the additional entries in the vector. The correspondence between problems on the two domains is hence only approximate, so our approach is one of preconditioning: using a fast Poisson solver on the rectangular region $D$ as a preconditioner for a conjugate gradient solver on $D\backslash\Omega$.

While undoubtedly more sophisticated schemes are possible, this approach to computing the biharmonic extension has a simple implementation, with the entire solution process described by the MATLAB code in Algorithm 1. We note that there is no need to form the 13-banded discrete biharmonic matrix at all with this iterative approach: the left hand side of (2.1) can simply be evaluated on demand.

The cost of applying the fast Poisson preconditioner is $O(N\log N)$ operations for a grid of $N=n\times n$ nodes in total. In the next section, we report for each of our numerical experiments the number of conjugate gradient iterations required for convergence using this preconditioner. Typically we found that the required number of iterations was a slowly growing function of the linear dimension $n$.

3 Numerical experiments

In this section we demonstrate and discuss the effectiveness of the biharmonic extension approach on five distinct examples, in two and three dimensions, including problems exhibiting symmetry and periodicity.

3.1 Example 1, two-dimensional peanut

For our first example, taken from [22], we consider the domain $D=[-\pi,\pi]\times[-\pi,\pi]$ with field values

[TABLE]

prescribed in the region

[TABLE]

This gives rise to the interface shown in black in Figure 2(a), with the field known on the inside. The goal is to smoothly extend the field to all of $D$. We also show in grey the extension given by formula (13) itself for reference. We will benchmark our computed extensions by comparing against this reference field in a small neighbourhood of the interface.

We test three different choices of boundary conditions on $D$, namely Dirichlet, Neumann and a mixture of Neumann on the sides and Dirichlet on the top and bottom. (Modifications to the code in Algorithm 1 for Neumann conditions are straightforward, requiring just a change to the array padding in the biharmonic stencil, and changing the eigenvalues and fast transforms from sine transforms to cosine transforms.) Figure 2(b)–(d) illustrates the biharmonic extension fields under these three choices of boundary conditions. The extension field shown in Figure 2(d) comes closest to representing the reference field in Figure 2(a), owing to the chosen boundary conditions matching the actual behaviour of formula (13). But very often in practice no such boundary information will be available, and the computational boundary is simply placed sufficiently far from the interface, as is the usual practice with level set methods. Hence the essential focus is on the behaviour of the extension in the neighbourhood of the interface and, to visual accuracy, the fields in Figure 2(b)–(d) do equally well in this regard.

Table 1 confirms this observation numerically, by tabulating the maximum difference between each extension and the reference field in a neighbourhood of four grid cells of the interface. As the number of grid divisions is doubled in each direction, the errors for all three choices of boundary condition are seen to reduce by approximately a factor of four, which is consistent with second order spatial accuracy.

Using the iterative solution strategy given in Algorithm 1, the number of conjugate gradient iterations required to compute the field on an $n\times n$ grid is seen to grow slightly faster than linearly with $n$, requiring for example 257 and 579 iterations for the mixed case on $512\times 512$ and $1024\times 1024$ grids respectively. All solutions were solved to iterative tolerance $10^{-6}$ in relative norm, which was sufficient for all tabulated errors to be unchanged with further iteration.

In fact, for this example the direct approach of simply forming and factorising the discrete biharmonic matrix proves to be more efficient than the iterative solver. In Table 2 we report the runtimes for both direct and iterative methods for the biharmonic extension, observing that the direct method is much faster for this problem. For example, with $n=1024$ the direct solver requires 12.4 seconds compared to 131.0 seconds for the iterative solver. For comparison purposes, we also report in Table 2 the errors and runtimes for computing constant, linear and quadratic extrapolations by solving (3), and its higher order variants, as described by Aslam [21]. As expected, these methods are observed to converge at first, second, and third order respectively to the reference field in a neighbourhood of the interface. Concerning runtime, the biharmonic extension computed iteratively is as fast, or faster, than all the extrapolation methods, while the biharmonic extension computed directly is a clear order of magnitude faster to compute than the other methods. We also report in Table 2 the condition numbers of the discrete biharmonic matrix, estimated using MATLAB’s condest function. These grow at the expected rate of 16 as the number of divisions is doubled, and are very similar to the condition numbers obtained from the squares of the eigenvalues of the discrete Laplacian on a rectangle, as in Algorithm 1.

The constant, linear and quadratic extrapolations are compared with the biharmonic extension (with Neumann conditions) in Figure 3, which reproduces the results presented in [22] Fig. 3. Although the biharmonic extension and linear extrapolation are of comparable accuracy in the vicinity of the interface, the biharmonic extension is free of the sharp corners that develop in the linear (and quadratic) extrapolation along the vertical centreline $x=0$, which are artifacts of the cusp in the interface at this location.

3.2 Example 2, three-dimensional peanut

For our second example, we consider an extension of Example 1 to three dimensions over the domain $D=[-\pi,\pi]^{3}$. Although the presentation to this point has concerned biharmonic extensions in two dimensions, all results generalise naturally to three. The radial basis function associated with the extension is now the linear RBF $\Phi(r)=r$ (which is a fundamental solution to the biharmonic equation in three dimensions) and (8) – (10) continue to hold unchanged. For this example we specify the field values

[TABLE]

prescribed in the region

[TABLE]

This gives rise to the peanut-shaped interface obtained by revolving the interface from Example 1 around the $x$ axis. Again, the field is known on the inside and the goal is to smoothly extend the field to all of $D$.

This time we consider periodic conditions on the biharmonic extension. In three dimensions, and especially with periodic conditions, the factorisation costs associated with the direct approach can be prohibitive owing to the amount of fill-in generated, so the iterative approach of Algorithm 1 is beneficial (the necessary changes are to use symmetric array padding and discrete Fourier transforms in all three dimensions.) Table 3 lists the runtimes, errors and estimated orders of convergence using this solution strategy. The number of conjugate gradient iterations required to compute the field on an $n\times n\times n$ grid is seen to grow only slightly faster than linearly with $n$, requiring for example 112 and 276 iterations on $128^{3}$ and $256^{3}$ grids respectively. By contrast, attempting to factorise the discrete biharmonic matrix quickly exhausted the available memory on the test machine. Statistics for the constant, linear and quadratic extrapolations are also reported in Table 3. The biharmonic extension is observed to be several times faster to compute than the comparable, linear extrapolation.

In Figure 4 we present contours for the $z\approx 0.6$ slice for this problem. All fields are accurate representations of the reference field in the vicinity of the interface. Away from the interface, the biharmonic extension is free from the steep gradients that are present especially in the linear and quadratic extrapolations.

3.3 Example 3, an annulus

For our third example we take the domain $D=[-2,2]\times[-2,2]$ with field values

[TABLE]

prescribed in the annulus

[TABLE]

The two interfaces defining the annulus, and the known field within it, are illustrated in Figure 5(a). Again, we also show in grey the contour lines for the reference field on all of $D$ for later comparison. Multiply-connected domains such as this pose no additional difficulty to the biharmonic extension. Indeed, the biharmonic extension is quite convenient in the sense that it extends the field to all of $D$ with a single solve (i.e. it extends both “inside” and “outside” at the same time). We note that although the reference field is singular at the origin, the computed extensions are well defined and smooth at all points in $D$.

Figure 5(b)–(c) illustrates the biharmonic extension fields subject to Dirichlet and Neumann boundary conditions. In the interior, the two extensions are identical, since the biharmonic stencils for interior nodes only reference other interior nodes or specified values within the annulus. On the outside, both extensions appear similar near the interface, with the Neumann conditions producing an extension more in keeping with the reference field away from the interface. Table 4 confirms that the error within four grid cells on the outside is in fact smaller for the Neumann extension than for the Dirichlet, but both extensions converge at the expected rate to the reference field as the grid spacing is halved. Once again, the number of conjugate gradient iterations grows slightly faster than linearly with $n$, up to a maximum of 522 iterations for $n=1024$ with Neumann boundary conditions.

For this problem, an alternative approach is to exploit the symmetry of the interface and the field about the $y$ axis. Applying Neumann conditions on the reduced problem over the rectangle $D=[0,2]\times[-2,2]$ produces the field shown in Figure 5(d): identical to that of Figure 5(c), but obtained using a grid with half as many nodes. The number of conjugate gradient iterations required for the $512\times 1024$ grid was also reduced compared to the full grid, down to 387 from 522.

The biharmonic extension (with Neumann conditions) is compared with constant, linear and quadratic extrapolations in Figure 6. Here we note that these extrapolations do inherit the singularity at the origin from the field equation (15), which the biharmonic extension does not. Nonetheless in the vicinity of the interface all computed fields exhibit their expected orders of convergence, as confirmed in Table 5. In this table we also report runtimes for all methods. For this two-dimensional problem, the direct method for the biharmonic extension is easily the fastest, at just 11.3 seconds, while the constant extrapolation and the iterative method for the biharmonic are of comparable runtime, at around 128 seconds. The linear and quadratic extrapolations require, as expected, two and three times longer to compute than this, respectively.

3.4 Example 4, periodic domain / channel

For our next example we take the domain $D=[0,3]\times[0,1]$ with periodic conditions on the top ($y=1$) and bottom ($y=0$). The field values

[TABLE]

are given on the left of the interface as shown in Figure 7(a). The interface itself is described by a sum of sine waves $\sum_{k=1}^{10}a_{k}\sin(2\pi ky+\omega_{k})$ with uniformly random amplitudes $a_{k}\in[-0.05,0.05]$ and uniformly random phase shifts $\omega_{k}\in[0,2\pi)$. This problem is designed to simulate an application where the stability of a travelling front to small sinusoidal perturbations is investigated. The periodic conditions are quite essential to the analysis, and hence the extension field must also respect this periodicity. The biharmonic extension naturally does so with the imposition of periodic conditions on the top and bottom, and we chose Neumann conditions for the sides.

Figure 7(b) shows the biharmonic extension for this problem, replicated over two additional periods for illustrative purposes (the problem is solved numerically over just one period). No special handling of nodes near the interface is required, as usual – the biharmonic extension naturally takes care of incorporating the required field information whenever a stencil crosses the interface. Close inspection of Figure 7(a) shows that some contour lines are cut off by the interface multiple times, alternately to the left and then to the right, but these and all other contours are smoothly extended in Figure 7(b). However, we do observe in Table 6 that the convergence rate of the biharmonic extension to the reference field is slower than in previous examples, estimated at only 1.5.

The lower convergence rate is due to the irregularity of the interface, resulting in biharmonic stencils near the peaks of the interface incorporating fewer nodes on the inside that carry known field values. In Figure 7(c), we see that the largest errors between the reference and extension fields occur in these regions near the peaks. Conversely, the errors are much lower in regions near the troughs of the interface, where the stencils reference many field values on the inside. If, instead, we use only a single period of a sine wave $a_{1}\sin(2\pi y+\omega_{1})$ as the shape of the interface (not pictured), we obtain the second set of results in Table 6. Here, convergence is restored to near second order, confirming that the irregularity of interface is the cause of the reduced convergence rate.

The preconditioned conjugate gradient iteration continues to be effective, although for the irregular interface on the finest grid with $N=1536\times 512$ nodes, the number of iterations required increased to more than $\sqrt{N}$, compared to less than $\sqrt{N}/2$ for the single period sine wave. This suggests that as grids are further refined, more advanced iterative solvers may be worth investigating. The direct solver is, naturally, unaffected by such concerns, and generates solutions more than an order of magnitude faster than the iterative solver on this two-dimensional problem.

Channel geometry

Alternatively, the problem may be considered in a channel, with Neumann conditions imposed on the top and bottom. For comparative purposes, we use the same interface, but now the prescribed field

[TABLE]

satisfies Neumann, rather than periodic, conditions at the top and bottom.

We present our results in Figure 8 and Table 7. There are no significant differences observed in the behaviours of the runtimes or the errors in Table 7 as compared to Table 6. From Figure 8, we observe that the field is smoothly extended across the interface, and is in keeping with the Neumann conditions prescribed on the boundaries. We note again that the error in the extended field is greatest near the peaks of the interface, due to the small number of neighbouring nodes that carry known field values. Indeed, for the coarsest grid, $192\times 64$, the biharmonic stencils at the tips of the steepest peaks reference just two values on the inside: $f_{i-1,j}$ and $f_{i-2,j}$, with all other neighbouring values being themselves constructed as part of the extension. Although this scenario leads to a local reduction in the accuracy, it causes no particular problems for the biharmonic extension in its implementation. In fact, it is worthwhile emphasising that even a single isolated node surrounded by the interface, or a collection of such isolated nodes, lead to well-posed biharmonic extensions. In the latter case the method would simply be building a clamped thin plate spline that interpolates the data at the isolated points, which seems as reasonable an extension as one could hope for under those circumstances.

3.5 Example 5, Stefan problem

For our final example, we implement the biharmonic extension in the context of a full level set method to solve a Stefan problem. The problem models solidification of a seed in an undercooled liquid, and is governed by the following dimensionless equations. Temperature in both solid ($s$) and liquid ($\ell$) phases obeys the heat equation,

[TABLE]

The normal velocity of the interface that separates the two phases is given by the jump condition

[TABLE]

The temperature on the interface is given by the Gibbs-Thomson relation,

[TABLE]

which models the effects of surface tension with coefficient $\sigma$, where $\kappa$ is the interface curvature. The model is subject to the following initial condition for temperature,

[TABLE]

where $\beta$ is the Stefan number. The temperature in the far-field is given by $T\rightarrow-1/\beta$ as $\|\mathbf{x}\|\rightarrow\infty$, which is imposed as Neumann boundary conditions on the computational domain $D$.

This problem was considered by Gibou et al. [14] who presented solutions based on a ghost-fluid level-set method previously proposed in [11]. Briefly, with the ghost-fluid method, the heat equation (18) is discretised using finite differences and, wherever computational stencils cross the interface, ghost values are constructed by extrapolation which implicitly satisfy the jump condition. Temperate gradients are calculated in each phase, and the normal derivative is then extrapolated in the normal direction, allowing a velocity field to be calculated at all grid points using (19). The level set method is then used to advance the moving front according to (2), with the level set function reinitialised to approximate a signed distance function after each step. Further details can be found in [11, 14, 20].

Our numerical scheme does not make use of ghost cells, and instead uses asymmetric stencils in the vicinity of the interface in the manner described by Chen et al. [10]. Hence we solve for only one temperature field at each node, that being the temperature of the phase present there. We compute gradients from these temperatures, again using asymmetric stencils in the vicinity of the interface, in both the $x$-$y$ and $\eta$-$\zeta$ (i.e. $45^{\circ}$) directions, as suggested by Chen et al. [10]. These gradient components, defined only at grid nodes, are then extended from one phase to the other using biharmonic extensions with Neumann boundary conditions. That is, the solid phase gradients are extended outwards using one biharmonic solve, and the liquid phase gradients are extended inwards using another biharmonic solve, allowing a velocity field to be calculated at all grid points using (19). From there the level set method proceeds as outlined above.

Figure 9 reproduces Fig. 24 of [14] using our approach. The Stefan number is held fixed at $\beta=2$, time runs to $t=0.4$, and the surface tension coefficient $\sigma$ is varied. Our solutions are shown for a $200\times 200$ grid. We used first order Euler time stepping schemes (implicit for the heat equation, explicit for the level set equation) with small fixed stepsize $\Delta t=0.0005$. We note that the Stefan problem with $\sigma=0$ is ill-posed, and the numerical scheme is in effect providing a form of regularisation, which allows Figure 9(a) to be obtained. Hence, with different numerical schemes it is not expected to obtain agreement for this figure, and indeed we find our Figure 9(a) is not identical to that published in [14], although both solutions are highly unstable. With $\sigma>0$ the problem is well-posed (although still unstable), and our numerical scheme incorporating the biharmonic extension produces results (Figure 9(b)-(c)) consistent with those published in [14].

These results serve to illustrate the utility of the biharmonic extension as part of a level set solution strategy. In our implementation, we utilised sparse direct factorisation in solving both the heat equation and the biharmonic equation. We found that factorising the discrete biharmonic matrix took two to three times as long as factorising the (scaled, shifted) discrete Laplacian matrix of the heat equation. This led to approximately 40% of the total runtime being taken up with biharmonic solves, 17% with heat equation solves, and the remaining runtime split between gradient calculations, level set reinitialisations, and other level set machinery.

4 Conclusion

In this paper we have presented a method for computing field extensions in level set methods by solving a biharmonic equation. This approach differs from other typical approaches in its fully implicit treatment of the interface. No explicit properties of the interface such as its location or the field value on the interface are required. This leads to a particularly simple implementation using either a sparse factorisation, or a matrix-free conjugate gradient solver with a fast Poisson preconditioner.

We demonstrated the biharmonic extension on a number of test problems that served to illustrate its effectiveness at producing smooth and accurate extensions near interfaces, whether extending inside or outside, on single or multiply-connected domains, and conforming to a variety of prescribed boundary conditions including symmetry and periodicity. As part of a level set solution strategy for moving interface problems, the method works well, as was demonstrated for an unstable Stefan problem. Further, we have shown there is no essential difference in applying this method in three dimensions, given that the biharmonic equation yields the natural choice of radial basis function $\Phi(r)=r$ in three dimensions, as it does with $\Phi(r)=r^{2}\log(r)$ in two.

Overall, we believe this approach will prove to be a very convenient method for incorporating velocity or other field extensions into level set methods. Developing further improvements to the efficiency of the iterative solution strategy, particularly important in three dimensions, seems an excellent avenue for future research.

Acknowledgements

TJM acknowledges the helpful discussions with Stephen Roberts from Australian National University. SWM acknowledges the support of the Australian Research Council via the Discovery Project DP140100933. All authors thank the anonymous referees for their suggestions that have helped improve this paper.