1 Introduction
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 levelset methods
[30], extrapolation procedures in space have been frequently used since the advent of the ghostfluid method [11], where constant extrapolations were originally used. Generalized ghostfluid methods were then designed, in part based on higherorder extrapolations for which Aslam introduced a partial differential equation (PDE) approach to perform linear and quadratic extrapolation
[2] and Gibou and Fedkiw introduced a cubic extrapolation in the same PDE framework [14]. It is natural in the levelset context to perform such extrapolations using PDE formulations for their solutions are based on HamiltonJacobi solvers that have been designed for other standard levelset equations, see e.g. [39]. 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 righthand 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 levelset applications including multiphase flow simulations [21, 13, 20, 28, 7, 17, 34], in the solution of PoissonBoltzmann [25, 19] and PoissonNerntzPlanck equations [23] for studying transport in ionic solutions, in heat and diffusion flow problems [14, 15, 4, 32], in the study of epitaxial growth and diblockcopolymer selfassembly used in the semiconductor industry [33, 31], in shape optimization [1, 40], surface reconstruction of biomolecules [25, 10] and in Stefantype problems [12, 6, 26, 35]. PDEbased extrapolation procedures have also been extended to adaptive Quad/Octree 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 [27] 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 highcurvature 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 Levelset Representation
The level set representation [30] defines the interface of a domain by , its interior and exterior by and , respectively, where is a Lipschitz continuous function called the levelset function. In this paper, the only geometrical quantity that is needed is the outward normal to the interface, , which can be computed as:
(1) 
using central differencing for and . In typical levelset simulations, the levelset function is reinitialized as a signed distance function [39]. We refer the interested reader to [36, 29] for a thorough presentation of the levelset method and [16] 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 [2]. 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:
(2) 
where is the Heaviside function. Then, the value of across the interface is found by solving the following two partial differential equations:
(3)  
(4) 
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 levelset region:
Similar to the method described in the previous section, we extrapolate the elements of in a constant fashion:
(5) 
before successively solving the following equations:
(6)  
(7) 
Note that now, the normal vector field
enters 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 [2] 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 secondorder 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 welldefined 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:
Applying an explicit firstorder accurate in time discretization to equations (2)(7) one obtains the following updating formulas:
(8) 
and
(9) 
Firstorder spatial derivatives are computed using secondorder upwind discretizations. For example, derivatives in the direction are approximated as:
where
Derivatives in the direction are approximated in a similar fashion.
Since for quadratic extrapolations it is sufficient to extend first and secondorder derivatives with firstorder accuracy only, the correction involving the minmod operator is dropped when solving the corresponding equations. Furthermore, since in the new method the approximation of secondorder 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 firstorder 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 steadystate 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 timedependent problems). Such an approach will be explored in future works.
2.4.1 Extension to adaptive Quad/Octree grids
The methodology introduced in this paper can be trivially extended to Quad/Octree data structures. Specifically, we sample data fields at nodes of Quad/Octree grids and use the xsecondorder accurate discretizations of [22] for regions where grids are nonuniform. 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 [2] performs well. The starshape 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 levelset 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 [2] 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 thirdorder rates of convergence for the linear and quadratic extrapolation, respectively.
For the highcurvature 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 wellresolved. 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 [2] produces large errors near kinks; those errors are significantly reduced with the Present approach. In particular, Figures 3cd and 4cd demonstrate that the secondorder (thirdorder) accuracy of the linear (quadratic) extrapolations are recovered with the proposed approach; the rates of convergence for the approach of [2] 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 highcurvature 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 levelset functions:
Figure 6 depicts those domains along with the octree grid refined near their boundaries.
Similar to the twodimensional 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 wellresolved 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.
5 Conclusion
We have presented a numerical method for extrapolating scalar quantities across the boundaries of irregular domains that may present highcurvature features or kinks. Linear and quadratic extrapolations procedures produce second and thirdorder 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 levelset methods. The numerical method we introduced is based on solving PDEs in pseudotime, 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.
Acknowledgment
This work was supported by ONR MURI N000141712676.
References
 [1] Grégoire Allaire, François Jouve, and AncaMaria Toader. Structural optimization using sensitivity analysis and a levelset mehtod. Journal of Computational Physics, 2003.
 [2] T Aslam. A partial differential equation approach to multidimensional extrapolation. J. Comput. Phys., 193:349–355, 2004.
 [3] 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.
 [4] 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.
 [5] A. Chacon and A. Vladimirsky. A parallel twoscale method for Eikonal equations. SIAM J. on Scientific Computing, 37(A156A180), 2015.
 [6] 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.
 [7] Charles Cleret de Langavant, Arthur Guittet, Maxime Theillard, Fernando TempranoColeto, and Frédéric Gibou. Levelset simulations of soluble surfactant driven flows. J. Comput. Phys., 348:271 – 297, 2017.
 [8] Miles Detrixhe and Frédéric Gibou. Hybrid Massively Parallel Fast Sweeping Method for static HamiltonJacobi Equations. In preparation, 2014.
 [9] 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.
 [10] Raphael Egan and Frédéric Gibou. Fast and scalable algorithms for constructing solventexcluded surfaces of large biomolecules. Journal of Computational Physics, 374:91 – 120, 2018.
 [11] Ronald P Fedkiw, Tariq Aslam, Barry Merriman, and Stanley Osher. A nonoscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method). J. Comput. Phys., 152(2):457 – 492, 1999.
 [12] 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.
 [13] 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.
 [14] 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.
 [15] Frederic Gibou, Ronald Fedkiw, L.T. Cheng, and Myngjoo Kang. A secondorder accurate symmetric discretization of the Poisson equation on irregular domains. J. Comput. Phys., 176:205–227, 2002.
 [16] Frederic Gibou, Ronald Fedkiw, and Stanley Osher. A review of levelset methods and some recent applications. J. Comp. Phys., 353:82–109, 2018.
 [17] Frederic Gibou and Chohong Min. Efficient symmetric positive definite secondorder accurate monolithic solver for fluid/solid interactions. J. Comp. Phys., 231:3246–3263, 2012.
 [18] 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.
 [19] Ásdís Helgadóttir and Frederic Gibou. A PoissonBoltzmann solver on irregular domains with Neumann or Robin boundary conditions on nongraded adaptive grid. J. Comput. Phys., 230:3830–3848, 2011.
 [20] Mathieu Lepilliez, Elena Roxana Popescu, Frederic Gibou, and Sébastien Tanguy. On twophase flow solvers in irregular domains with contact line. J. Comput. Phys., 321:1217 – 1251, 2016.
 [21] 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.
 [22] Chohong Min and Frederic Gibou. A second order accurate level set method on nongraded adaptive Cartesian grids. J. Comput. Phys., 225(1):300–321, 2007.
 [23] 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.
 [24] M. Mirzadeh, A. Guittet, C. Burstedde, and F. Gibou. Parallel levelset methods on adaptive treebased grids. J. Comp. Phys. (submitted), 2015.
 [25] M. Mirzadeh, M. Theillard, A. Helgadottir, D. Boy, and F. Gibou. An adaptive, finite difference solver for the nonlinear PoissonBoltzmann equation with applications to biomolecular computations. Communications in Computational Physics, 13(1):150–173, 2012.
 [26] 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.
 [27] 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.
 [28] 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.
 [29] S. Osher and R. Fedkiw. Level Set Methods and Dynamic Implicit Surfaces. SpringerVerlag, 2002. New York, NY.
 [30] Stanley Osher and James A. Sethian. Fronts propagating with curvature dependent speed: algorithms based on HamiltonJacobi formulations. J. Comput. Phys., 79(1):12–49, 1988.
 [31] Gaddiel Y. Ouaknin, Nabil Laachi, Kris Delaney, Glenn H. Fredrickson, and Frederic Gibou. Levelset strategy for inverse DSAlithography. Under review, 2017.
 [32] Joseph Papac, Frederic Gibou, and Christian Ratsch. Efficient symmetric discretization for the Poisson, heat and Stefantype problems with Robin boundary conditions. J. Comput. Phys., 229:875–889, 2010.
 [33] Joseph Papac, Asdis Helgadottir, Christian Ratsch, and Frederic Gibou. A level set approach for diffusion and Stefantype problems with Robin boundary conditions on Quadtree/Octree adaptive Cartesian grids. J. Comput. Phys., 233:241–261, 2013.
 [34] A. RobinsonMosher, T. Shinar, J. Gretarsson, J. Su, and R. Fedkiw. Twoway coupling of fluids to rigid and deformable solids and shells. ACM Trans. Graph., 27(46), 2008.
 [35] 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.

[36]
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.
 [37] J.A. Sethian and A. Vladimirsky. Fast methods for the Eikonal and related HamiltonJacobi equations on unstructured meshes. Proc. Natl. Acad. Sci., 97/11:5699–5703, 2000.
 [38] J.A. Sethian and A. Vladimirsky. Ordered upwind methods for static hamiltonjacobi equations. Proc. Natl. Acad. Sci, 98/20:11069–11074, 2001.
 [39] M Sussman, P Smereka, and S Osher. A level set approach for computing solutions to incompressible twophase flow. J. Comp. Phys., 114:146–159, 1994.
 [40] Maxime Theillard, Landry Fokoua Djodom, JeanLéopold Vié, and Frédéric Gibou. A secondorder 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.
 [41] Hongkai Zhao. A fast sweeping method for eikonal equations. Mathematics of Computation, 74:603–627, 2004.
Comments
There are no comments yet.