Extrapolation procedures are ubiquitous in scientific computing and generally allow one to estimate a valid value of a quantity at points where data is not given; either in space or in time. In the context of level-set methods, extrapolation procedures in space have been frequently used since the advent of the ghost-fluid method 
, where constant extrapolations were originally used. Generalized ghost-fluid methods were then designed, in part based on higher-order extrapolations for which Aslam introduced a partial differential equation (PDE) approach to perform linear and quadratic extrapolation and Gibou and Fedkiw introduced a cubic extrapolation in the same PDE framework . It is natural in the level-set context to perform such extrapolations using PDE formulations for their solutions are based on Hamilton-Jacobi solvers that have been designed for other standard level-set equations, see e.g. . A typical situation that needs extrapolation is that of an implicit treatment of a field in a free boundary problem. In this case, a valid value of the field at time needs to be known when assembling the right-hand side of the linear system of equations at time . Since the interface at the new time step has swept grid points that are outside the domain at the previous time step, valid values of the field at time are needed in the domain at time , which requires an extrapolation procedure.
Typical use of extrapolation methods can be found in a multitude of level-set applications including multiphase flow simulations [21, 13, 20, 28, 7, 17, 34], in the solution of Poisson-Boltzmann [25, 19] and Poisson-Nerntz-Planck equations  for studying transport in ionic solutions, in heat and diffusion flow problems [14, 15, 4, 32], in the study of epitaxial growth and diblock-copolymer self-assembly used in the semi-conductor industry [33, 31], in shape optimization [1, 40], surface reconstruction of biomolecules [25, 10] and in Stefan-type problems [12, 6, 26, 35]. PDE-based extrapolation procedures have also been extended to adaptive Quad-/Oc-tree grids and parallel architectures [22, 24, 18]. In addition, fast methods have been introduced for computationally efficient extrapolation procedures using the Fast Marching method, including parallel implementations [37, 38, 5] or the Fast Sweeping method [41, 3], including efficient parallel algorithms on adaptive grids [9, 8]. We also refer the interested reader to  for another implicit approach to extrapolation based on solving the biharmonic equation.
However, those methods behave poorly in the case where the free boundary presents high-curvature features or kinks, which occur naturally in front propagation simulations, e.g. in the case where thin dendrites develop or when fronts merge. We introduce a method that solves that problem. We present the method in section 2 and numerical examples in sections 3 and 4 that illustrate its benefits and comment on its efficiency. Section 5 draws some conclusions.
2 Numerical Method
2.1 Level-set Representation
The level set representation  defines the interface of a domain by , its interior and exterior by and , respectively, where is a Lipschitz continuous function called the level-set function. In this paper, the only geometrical quantity that is needed is the outward normal to the interface, , which can be computed as:
using central differencing for and . In typical level-set simulations, the level-set function is reinitialized as a signed distance function . We refer the interested reader to [36, 29] for a thorough presentation of the level-set method and  for a recent review.
2.2 Standard multidimensional extrapolation
High order extrapolations in the normal direction are traditionally performed in a series of steps, as proposed by Aslam in . For example, suppose that we seek to extrapolate a scalar field from the region where to the region where . In the case of a quadratic extrapolation, we first compute in the region and extrapolate it across the interface in a constant fashion, that is, such that its normal derivative is zero in the region , by solving the following partial differential equation:
where is the Heaviside function. Then, the value of across the interface is found by solving the following two partial differential equations:
defining in such a way that its normal derivative is equal to the previously extrapolated and then defining in such a way that its normal derivative is equal to the previously extrapolated . These PDEs are solved in fictitious time for a few iterations (typically 15) since we only seek to extrapolate the values of in a narrow band of a few grid cells around the interface.
2.3 Present multidimensional extrapolation
Instead of calculating the normal derivatives in the negative region before extrapolating them, we instead compute the derivatives in the Cartesian directions, extrapolate them and then construct the normal derivatives. Specifically, consider the following quantities, that are computed in the negative level-set region:
Similar to the method described in the previous section, we extrapolate the elements of in a constant fashion:
before successively solving the following equations:
Note that now, the normal vector fieldenters the equations merely as some sort of weighting factor. Thus, as long as field is sufficiently smooth this approach to multidimensional extrapolation is expected to produce accurate results even when the normal vector field is not smooth (as is the case of domains with sharp features).
2.4 Implementation details
In this work we demonstrate the proposed method on uniform Cartesian grids and our implementation follows very closely the one from  with just few differences. Consider a two dimensional computational grid with nodes defined as:
where denotes the computational domain, and are number of grid nodes in the Cartesian directions. Standard second-order accurate central difference formulas are used for calculating the normal vector field (in the entire domain) and derivatives (first and second) of in the negative region. Normal derivatives of are computed as:
Since the first and second order derivatives of are not well-defined at all grid points where we replace the Heaviside function in equations (3), (6) and in equations (2), (5) with discrete fields and , respectively, where:
First-order spatial derivatives are computed using second-order upwind discretizations. For example, derivatives in the -direction are approximated as:
Derivatives in the -direction are approximated in a similar fashion.
Since for quadratic extrapolations it is sufficient to extend first- and second-order derivatives with first-order accuracy only, the correction involving the minmod operator is dropped when solving the corresponding equations. Furthermore, since in the new method the approximation of second-order derivatives in all Cartesian directions are already available during solving the PDE for , this correction can be computed only once during first iteration and reused in subsequent iterations, reducing the cost of each iteration by approximately 2 times. Specifically, the total count of arithmetic operations to compute using the Standard method is approximately 22 in two spatial dimensions and 32 in three spatial dimensions, while for the Present method the total count is 10 and 14, correspondingly. Thus, if we denote as the cost of solving a single advection equation using first-order accurate approximations of derivatives, then the total cost of performing quadratic extrapolation using the Standard method is approximately in two and three spatial dimensions, while the total cost of performing quadratic extrapolation using the Present method is in two spatial dimensions and in three spatial dimensions.
Remark. Since in the Present approach there is no need to recalculate second derivatives and apply the nonlinear minmod operator at every iteration, it is possible to obtain the steady-state solution of the advection equations in an implicit fashion. This could be very beneficial in cases when a good guess for the extended field is available (for example, solutions from preceding time instants in time-dependent problems). Such an approach will be explored in future works.
2.4.1 Extension to adaptive Quad-/Oc-tree grids
The methodology introduced in this paper can be trivially extended to Quad-/Oc-tree data structures. Specifically, we sample data fields at nodes of Quad-/Oc-tree grids and use the xsecond-order accurate discretizations of  for regions where grids are non-uniform. A band of uniform grids are usually imposed near the interface in practical free boundary applications (see Fig. 1). In this case the extrapolation within some neighborhood around the interface (where it is primarily required) is as accurate as for uniform grids, however, the extrapolation procedure is much faster on adaptive grids for their significant reduction in the total number of grid points.
3 Numerical Results in Two Spatial Dimensions
We consider four physical domains: a disk, a star shape, a union of two disks and an intersection of two disks (see figure 2). The disk is a smooth interface for which the Standard approach of  performs well. The star-shape domain is an example where regions of high curvature are present (crest and trough of the wavy shape). The union/intersection of two disks are examples where kinks occur and illustrate the case of typical free boundary simulations where changes in topology occur. The definition of those domains are given by the level-set functions:
In each case we consider a computational domain . We define the function inside every domain and extrapolate it in the outside region. Then the maximum difference between the exact values of and extrapolated ones is computed within a band of thickness in the outside regions. Figures 3 and 4 summarize the convergence behavior of the Standard approach of  and the Present one. Figure 5 demonstrates the error distribution for both methods in the case of the quadratic extrapolation on a grid.
In case of the smooth domain (disk), both approaches produce almost indistinguishable results attaining second- and third-order rates of convergence for the linear and quadratic extrapolation, respectively.
For the high-curvature domain (star), both methods still reach optimal orders of convergence; however the Standard approach demonstrates the optimal order of convergence only at relatively high grid resolutions when all geometric features are well-resolved. Moreover, for a given grid resolution the Present approach produces results that are more than one order of magnitude more accurate in the case of the linear extrapolation and almost three orders of magnitude more accurate in the case of the quadratic extrapolation compared to the Standard approach. Figure 5b shows that the Standard approach produces very large errors near regions with the highest curvature, while the error in the case of the Present approach is much smaller and exhibits very little variation throughout all regions around the interface.
The results are even more significantly improved with the proposed approach in the case of interfaces with kinks (union) and (intersection). Figures 5c and 5d show that the Standard method of  produces large errors near kinks; those errors are significantly reduced with the Present approach. In particular, Figures 3c-d and 4c-d demonstrate that the second-order (third-order) accuracy of the linear (quadratic) extrapolations are recovered with the proposed approach; the rates of convergence for the approach of  are close to first order, which corresponds to the constant extrapolation, due to the fact that errors near kinks do not decrease despite grid refinement.
4 Numerical Results in Three Spatial Dimensions
We consider three different domains, , and , that present high-curvature features or kinks in three spatial dimensions. In addition, we consider a smooth spherical domain with center and radius . The definition of those domains are given by the level-set functions:
Figure 6 depicts those domains along with the octree grid refined near their boundaries.
Similar to the two-dimensional examples, we consider a computational domain . We extrapolate the function from the inside to the outside for every domain and compute the difference between the exact values of and the extrapolated ones within a band of thickness in the outside region.
Conclusions similar to the two dimensional case can be drawn from the results in figures 7 and 8. Specifically, for a smooth and well-resolved domain (sphere) both approaches produce almost indistinguishable results with optimal order of convergence (second and third for the linear and quadratic extrapolations, respectively). When the interface curvature is high (, star) the Present approach produce extrapolated fields that are several orders of magnitude more accurate than for the Standard approach. For geometries with sharp features (union) and (intersection) only the Present approach demonstrates optimal orders of convergence, while for the Standard approach the rate of convergence is stuck to 1.
We have presented a numerical method for extrapolating scalar quantities across the boundaries of irregular domains that may present high-curvature features or kinks. Linear and quadratic extrapolations procedures produce second- and third-order accurate results in the norm, respectively and do so regardless of the irregularity of the boundaries, i.e. boundaries with kinks can readily be considered. These procedures are effective in both two and three spatial dimensions and can be implemented on quadtree and octree Cartesian grids. We have shown through numerical examples that errors associated with extrapolations can be reduced by several orders magnitude in some cases, compared with the current Standard approach used in level-set methods. The numerical method we introduced is based on solving PDEs in pseudo-time, but we note that static solutions based on an implicit approach like Fast Marching or Fast Sweeping could be obtained and we expect the results to follow the same general behavior as that presented in the current manuscript.
This work was supported by ONR MURI N00014-17-1-2676.
-  Grégoire Allaire, François Jouve, and Anca-Maria Toader. Structural optimization using sensitivity analysis and a level-set mehtod. Journal of Computational Physics, 2003.
-  T Aslam. A partial differential equation approach to multidimensional extrapolation. J. Comput. Phys., 193:349–355, 2004.
-  Tariq Aslam, songting Luo, and Hongkai Zha. Static pde approach for multidimensional extrapolation using fast sweeping methods. SIAM J. Sci. Comput., 36:A2907–A2928, 2014.
-  Daniil Bochkov and Frederic Gibou. Solving the poisson equation with Robin boundary conditions on piecewise smooth irregular boundaries. J. Comp. Phys., 376:1156–1198, 2019.
-  A. Chacon and A. Vladimirsky. A parallel two-scale method for Eikonal equations. SIAM J. on Scientific Computing, 37(A156-A180), 2015.
-  Han Chen, Chohong Min, and Frederic Gibou. A numerical scheme for the Stefan problem on adaptive Cartesian grids with supralinear convergence rate. J. Comput. Phys., 228(16):5803–5818, 2009.
-  Charles Cleret de Langavant, Arthur Guittet, Maxime Theillard, Fernando Temprano-Coleto, and Frédéric Gibou. Level-set simulations of soluble surfactant driven flows. J. Comput. Phys., 348:271 – 297, 2017.
-  Miles Detrixhe and Frédéric Gibou. Hybrid Massively Parallel Fast Sweeping Method for static Hamilton-Jacobi Equations. In preparation, 2014.
-  Miles Detrixhe, Frédéric Gibou, and Chohong Min. A parallel fast sweeping method for the eikonal equation. J. Comput. Phys., 237:46 – 55, 2013.
-  Raphael Egan and Frédéric Gibou. Fast and scalable algorithms for constructing solvent-excluded surfaces of large biomolecules. Journal of Computational Physics, 374:91 – 120, 2018.
-  Ronald P Fedkiw, Tariq Aslam, Barry Merriman, and Stanley Osher. A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method). J. Comput. Phys., 152(2):457 – 492, 1999.
-  F. Gibou, R. Fedkiw, R. Caflisch, and S. Osher. A level set approach for the numerical simulation of dendritic growth. J. Sci. Comput., 19:183–199, 2003.
-  Frederic Gibou, Liguo Chen, Duc Nguyen, and Sanjoy Banerjee. A level set based sharp interface method for the multiphase incompressible navier–stokes equations with phase change. J. Comput. Phys., 222(2):536–555, March 2007.
-  Frédéric Gibou and Ronald Fedkiw. A fourth order accurate discretization for the Laplace and heat equations on arbitrary domains, with applications to the Stefan problem. J. Comput. Phys., 202(2):577 – 601, 2005.
-  Frederic Gibou, Ronald Fedkiw, L.-T. Cheng, and Myngjoo Kang. A second-order accurate symmetric discretization of the Poisson equation on irregular domains. J. Comput. Phys., 176:205–227, 2002.
-  Frederic Gibou, Ronald Fedkiw, and Stanley Osher. A review of level-set methods and some recent applications. J. Comp. Phys., 353:82–109, 2018.
-  Frederic Gibou and Chohong Min. Efficient symmetric positive definite second-order accurate monolithic solver for fluid/solid interactions. J. Comp. Phys., 231:3246–3263, 2012.
-  Frederic Gibou, Chohong Min, and Ronald Fedkiw. High resolution sharp computational methods for elliptic and parabolic problems in complex geometries. J. Sci. Comput., 54:369–413, 2013.
-  Ásdís Helgadóttir and Frederic Gibou. A Poisson-Boltzmann solver on irregular domains with Neumann or Robin boundary conditions on non-graded adaptive grid. J. Comput. Phys., 230:3830–3848, 2011.
-  Mathieu Lepilliez, Elena Roxana Popescu, Frederic Gibou, and Sébastien Tanguy. On two-phase flow solvers in irregular domains with contact line. J. Comput. Phys., 321:1217 – 1251, 2016.
-  Frank Losasso, Frederic Gibou, and Ron Fedkiw. Simulating water and smoke with an octree data structure. ACM Trans. Graph. (SIGGRAPH Proc.), pages 457–462, 2004.
-  Chohong Min and Frederic Gibou. A second order accurate level set method on non-graded adaptive Cartesian grids. J. Comput. Phys., 225(1):300–321, 2007.
-  M. Mirzadeh and F. Gibou. A conservative discretization of the Poisson–Nernst–Planck equations on adaptive cartesian grids. J. Comput. Phys., 274:633–653, 2014.
-  M. Mirzadeh, A. Guittet, C. Burstedde, and F. Gibou. Parallel level-set methods on adaptive tree-based grids. J. Comp. Phys. (submitted), 2015.
-  M. Mirzadeh, M. Theillard, A. Helgadottir, D. Boy, and F. Gibou. An adaptive, finite difference solver for the nonlinear Poisson-Boltzmann equation with applications to biomolecular computations. Communications in Computational Physics, 13(1):150–173, 2012.
-  Pouria Mistani, Arthur Guittet, Daniil Bochkov, Joshua Schneider, Dionisios Margetis, Christian Ratsch, and Frederic Gibou. The island dynamics model on parallel quadtree grids. In preparation, 2017.
-  Timothy J Moroney, Dylan R Lusmore, Scott W McCue, and DL Sean McElwain. Extending fields in a level set method by solving a biharmonic equation. J. Comput. Phys., 343:170–185, 2017.
-  Duc Nguyen, Frederic Gibou, and Ronald Fedkiw. A Fully Conservative Ghost Fluid Method and Stiff Detonation Waves. In 12th Int. Detonation Symposium, San Diego, CA, 2002.
-  S. Osher and R. Fedkiw. Level Set Methods and Dynamic Implicit Surfaces. Springer-Verlag, 2002. New York, NY.
-  Stanley Osher and James A. Sethian. Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations. J. Comput. Phys., 79(1):12–49, 1988.
-  Gaddiel Y. Ouaknin, Nabil Laachi, Kris Delaney, Glenn H. Fredrickson, and Frederic Gibou. Level-set strategy for inverse DSA-lithography. Under review, 2017.
-  Joseph Papac, Frederic Gibou, and Christian Ratsch. Efficient symmetric discretization for the Poisson, heat and Stefan-type problems with Robin boundary conditions. J. Comput. Phys., 229:875–889, 2010.
-  Joseph Papac, Asdis Helgadottir, Christian Ratsch, and Frederic Gibou. A level set approach for diffusion and Stefan-type problems with Robin boundary conditions on Quadtree/Octree adaptive Cartesian grids. J. Comput. Phys., 233:241–261, 2013.
-  A. Robinson-Mosher, T. Shinar, J. Gretarsson, J. Su, and R. Fedkiw. Two-way coupling of fluids to rigid and deformable solids and shells. ACM Trans. Graph., 27(46), 2008.
-  Chris H. Rycroft and Frédéric Gibou. Simulations of a stretching bar using a plasticity model from the shear transformation zone theory. J. Comput. Phys., 231(5):2155 – 2179, 2012.
J. A. Sethian.
Level set methods, volume 3 of Cambridge Monographs on
Applied and Computational Mathematics.
Cambridge University Press, Cambridge, 1996.
Evolving interfaces in geometry, fluid mechanics, computer vision, and materials science.
-  J.A. Sethian and A. Vladimirsky. Fast methods for the Eikonal and related Hamilton-Jacobi equations on unstructured meshes. Proc. Natl. Acad. Sci., 97/11:5699–5703, 2000.
-  J.A. Sethian and A. Vladimirsky. Ordered upwind methods for static hamilton-jacobi equations. Proc. Natl. Acad. Sci, 98/20:11069–11074, 2001.
-  M Sussman, P Smereka, and S Osher. A level set approach for computing solutions to incompressible two-phase flow. J. Comp. Phys., 114:146–159, 1994.
-  Maxime Theillard, Landry Fokoua Djodom, Jean-Léopold Vié, and Frédéric Gibou. A second-order sharp numerical method for solving the linear elasticity equations on irregular domains and adaptive grids – application to shape optimization. J. Comput. Phys., 233:430 – 448, 2013.
-  Hongkai Zhao. A fast sweeping method for eikonal equations. Mathematics of Computation, 74:603–627, 2004.