Computation of incompressible bubble dynamics with a stabilized finite element level set method

This paper presents a stabilized finite element method for the three dimensional computation of incompressible bubble dynamics using a level set method. The interface between the two phases is resolved using the level set approach developed by Sethian [Level Set Methods and Fast Marching Methods, Cambridge University Press, 1999], Sussman et al. [J. Comput. Phys. 114 (1994) 146], and Sussman et al. [J. Comput. Phys. 148 (1999) 81–124]. In this approach the interface is represented as a zero level set of a smooth function. The streamline-upwind/Petrov–Galerkin method was used to discretize the governing flow and level set equations. The continuum surface force (CSF) model proposed by Brackbill et al. [J. Comput. Phys. 100 (1992) 335–354] was applied in order to account for surface tension effects. To restrict the interface from moving while re-distancing, an improved re-distancing scheme proposed in the finite difference context [J. Comput. Phys. 148 (1999) 81–124] is adapted for finite element discretization. This enables us to accurately compute the flows with large density and viscosity differences, as well as surface tension. The capability of the resultant algorithm is demonstrated with two and three dimensional numerical examples of a single bubble rising through a quiescent liquid, and two bubble coalescence.     

A New Level-Set Model for the Representation of Non-Smooth Geometries

In Starinshak et al. (J Comput Phys 262(1):1–16, 2014), we proposed a new level-set model for representing multimaterial flows in multiple space dimensions. Rather than associating each level-set function with the boundary of a material, the new model associates each level-set function with a pair of materials and the interface that separates them. In this paper, we extend the model to represent geometries with non-smooth boundaries. The model uses multiple level-set functions to describe the shape boundary, typically with one level-set function per smooth boundary segment. Sign information is collected from all level-set functions and a voting algorithm is used to determine the interior/exterior of the geometric shape. The model is well suited for representing boundaries with singularities; it offers significant improvement over standard level-set approaches, both in shape preservation and area conservation; and it eliminates the need for costly redistancing of the level-set function. Numerical examples illustrate the superior performance of the proposed model.     

Spatially adaptive techniques for level set methods and incompressible flow

Since the seminal work of [Sussman, M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114:146–59] on coupling the level set method of [Osher S, Sethian J. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations. J Comput Phys 1988;79:12–49] to the equations for two-phase incompressible flow, there has been a great deal of interest in this area. That work demonstrated the most powerful aspects of the level set method, i.e. automatic handling of topological changes such as merging and pinching, as well as robust geometric information such as normals and curvature. Interestingly, this work also demonstrated the largest weakness of the level set method, i.e. mass or information loss characteristic of most Eulerian capturing techniques. In fact, [Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114:146–59] introduced a partial differential equation for battling this weakness, without which their work would not have been possible. In this paper, we discuss both historical and most recent works focused on improving the computational accuracy of the level set method focusing in part on applications related to incompressible flow due to both of its popularity and stringent accuracy requirements. Thus, we discuss higher order accurate numerical methods such as Hamilton–Jacobi WENO [Jiang G-S, Peng D. Weighted ENO schemes for Hamilton–Jacobi equations. SIAM J Sci Comput 2000;21:2126–43], methods for maintaining a signed distance function, hybrid methods such as the particle level set method [Enright D, Fedkiw R, Ferziger J, Mitchell I. A hybrid particle level set method for improved interface capturing. J Comput Phys 2002;183:83–116] and the coupled level set volume of fluid method [Sussman M, Puckett EG. A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows. J Comput Phys 2000;162:301–37], and adaptive gridding techniques such as the octree approach to free surface flows proposed in [Losasso F, Gibou F, Fedkiw R. Simulating water and smoke with an octree data structure, ACM Trans Graph (SIGGRAPH Proc) 2004;23:457–62].     

An analysis of and a comparison between the discontinuous Galerkin and the spectral finite volume methods

The discontinuous Galerkin method has been developed and applied extensively to solve hyperbolic conservation laws in recent years. More recently Wang et al. developed a class of discontinuous Petrov–Galerkin method, termed spectral (finite) volume method [J. Comput. Phys. 78 (2002) 210; J. Comput. Phys. 179 (2002) 665; J. Sci. Comput. 20 (2004) 137]. In this paper we perform a Fourier type analysis on both methods when solving linear one-dimensional conservation laws. A comparison between the two methods is given in terms of accuracy, stability, and convergence. Numerical experiments are performed to validate this analysis and comparison.     

Hybrid Finite Difference Weighted Essentially Non-oscillatory Schemes for the Compressible Ideal Magnetohydrodynamics Equation

In this paper, we present hybrid weighted essentially non-oscillatory (WENO) schemes with several discontinuity detectors for solving the compressible ideal magnetohydrodynamics (MHD) equation. Li and Qiu (J Comput Phys 229:8105–8129, 2010) examined effectiveness and efficiency of several different troubled-cell indicators in hybrid WENO methods for Euler gasdynamics. Later, Li et al. (J Sci Comput 51:527–559, 2012) extended the hybrid methods for solving the shallow water equations with four better indicators. Hybrid WENO schemes reduce the computational costs, maintain non-oscillatory properties and keep sharp transitions for problems. The numerical results of hybrid WENO-JS/WENO-M schemes are presented to compare the ability of several troubled-cell indicators with a variety of test problems. The focus of this paper, we propose optimal and reliable indicators for performance comparison of hybrid method using troubled-cell indicators for efficient numerical method of ideal MHD equations. We propose a modified ATV indicator that uses a second derivative. It is advantageous for differential discontinuity detection such as jump discontinuity and kink. A detailed numerical study of one-dimensional and two-dimensional cases is conducted to address efficiency (CPU time reduction and more accurate numerical solution) and non-oscillatory property problems. We demonstrate that the hybrid WENO-M scheme preserves the advantages of WENO-M and the ratio of computational costs of hybrid WENO-M and hybrid WENO-JS is smaller than that of WENO-M and WENO-JS.     

Weights Design For Maximal Order WENO Schemes

Weighted essentially non-oscillatory (WENO) finite difference schemes, developed by Liu et al. (Comput Phys 115(1):200–212, 1994) and improved by Jiang and Shu (Comput Phys 126(1):202–228, 1996), are one of the most popular methods to approximate the solutions of hyperbolic equations. But these schemes fail to provide maximal order accuracy near smooth extrema, where the first derivative of the solution becomes zero. Some authors have addressed this problem with different weight designs. In this paper we focus on the weights proposed by Yamaleev and Carpenter (J Comput Phys 228:4248–4272, 2009). They propose new weights to provide faster weight convergence than those presented in Borges et al. (J Comput Phys 227:3191–3211, 2008) and deduce some constraints on the weights parameters to guarantee that the WENO scheme has maximal order for sufficiently smooth solutions with an arbitrary number of vanishing derivatives. We analyze the scheme with the weights proposed in Yamaleev and Carpenter (J Comput Phys 228:4248–4272, 2009) and prove that near discontinuities it achieves worse orders than classical WENO schemes. In order to solve these accuracy problems, we define new weights, based on those proposed in Yamaleev and Carpenter (J Comput Phys 228:4248–4272, 2009), and get some constraints on the weights parameters to guarantee maximal order accuracy for the resulting schemes.     

Comparison and validation of compressible flow simulations of laser-induced cavitation bubbles

The numerical simulation of compressible two-phase fluid flows exhibits severe difficulties, in particular, when strong strong variations in the material parameters and high interface velocities are present at the phase boundary. Although several models and discretizations have been developed in the past, a thorough quantitative validation by experimental data and a detailed comparison of numerical schemes are hardly available.Here, two different discretizations are investigated, namely, a non-conservative approach proposed by Saurel and Abgrall [SIAM J Sci Comput 1999;21:1115] and the real ghost fluid method developed by Tang et al. [SIAM J Sci Comput 2006;28:278]. The validation is performed for the case of laser-induced cavitation bubbles collapsing in an infinite medium. For the computations, initial data are deduced implicitly from the experimental data. In particular, the influence of numerical phase transition caused by smearing of the phase boundary is investigated.     

A New Mapped Weighted Essentially Non-oscillatory Scheme

The weighted essentially non-oscillatory (WENO) methods are a popular high-order spatial discretization for hyperbolic partial differential equations. Recently Henrick et al. (J. Comput. Phys. 207:542–567, 2005) noted that the fifth-order WENO method by Jiang and Shu (J. Comput. Phys. 126:202–228, 1996) is only third-order accurate near critical points of the smooth regions in general. Using a simple mapping function to the original weights in Jiang and Shu (J. Comput. Phys. 126:202–228, 1996), Henrick et al. developed a mapped WENO method to achieve the optimal order of accuracy near critical points. In this paper we study the mapped WENO scheme and find that, when it is used for solving the problems with discontinuities, the mapping function in Henrick et al. (J. Comput. Phys. 207:542–567, 2005) may amplify the effect from the non-smooth stencils and thus cause a potential loss of accuracy near discontinuities. This effect may be difficult to be observed for the fifth-order WENO method unless a long time simulation is desired. However, if the mapping function is applied to seventh-order WENO methods (Balsara and Shu in J. Comput. Phys. 160:405–452, 2000), the error can increase much faster so that it can be observed with a moderate output time. In this paper a new mapping function is proposed to overcome this potential loss of accuracy.     

A Conservative Method for the Simulation of?the?Isothermal Euler System with the van-der-Waals Equation of State

In this article, we are interested in the simulation of phase transition in compressible flows, with the isothermal Euler system, closed by the van-der-Waals model. We formulate the problem as an hyperbolic system, with a source term located at the interface between liquid and vapour. The numerical scheme is based on (Abgrall and Saurel, J. Comput. Phys. 186(2):361?C396, 2003; Le Métayer et al., J. Comput. Phys. 205(2):567?C610, 2005). Compared with previous discretizations of the van-der-Waals system, the novelty of this algorithm is that it is fully conservative. Its Godunov-type formulation allows an easy implementation on multi-dimensional unstructured meshes.     

A new 3D parallel SPH scheme for free surface flows

We propose a new robust and accurate SPH scheme, able to track correctly complex three-dimensional non-hydrostatic free surface flows and, even more important, also able to compute an accurate and little oscillatory pressure field. It uses the explicit third order TVD Runge-Kutta scheme in time, following Shu and Osher [Shu C-W, Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J Comput Phys 1988;89:439-71], together with the new key idea of introducing a monotone upwind flux for the density equation, thus removing any artificial viscosity term. For the discretization of the velocity equation, the non-diffusive central flux has been used. A new flexible approach to impose the boundary conditions at solid walls is also proposed. It can handle any moving rigid body with arbitrarily irregular geometry. It does neither produce oscillations in the fluid pressure in proximity of the interfaces, nor does it have a restrictive impact on the stability condition of the explicit time stepping method, unlike the repellent boundary forces of Monaghan [Monaghan JJ. Simulating free surface flows with SPH. J Comput Phys 1994;110:399-406]. To asses the accuracy of the new SPH scheme, a 3D mesh-convergence study is performed for the strongly deforming free surface in a 3D dam-break and impact-wave test problem providing very good results.Moreover, the parallelization of the new 3D SPH scheme has been carried out using the message passing interface (MPI) standard, together with a dynamic load balancing strategy to improve the computational efficiency of the scheme. Thus, simulations involving millions of particles can be run on modern massively parallel supercomputers, obtaining a very good performance, as confirmed by a speed-up analysis. The 3D applications consist of environmental flow problems, such as dam-break flows and impact flows against a wall. The numerical solutions obtained with our new 3D SPH code have been compared with either experimental results or with other numerical reference solutions, obtaining in all cases a very satisfactory agreement.     

Symmetric and symplectic ERKN methods for oscillatory Hamiltonian systems

The ERKN methods proposed by H. Yang et al. [Comput. Phys. Comm. 180 (2009) 1777] are an important improvement of J.M. Franco?s ARKN methods for perturbed oscillators [J.M. Franco, Comput. Phys. Comm. 147 (2002) 770]. This paper focuses on the symmetry and symplecticity conditions for ERKN methods solving oscillatory Hamiltonian systems. Two examples of symmetric and symplectic ERKN (SSERKN) methods of orders two and four respectively are constructed. The phase and stability properties of the new methods are analyzed. The results of numerical experiments show the robustness and competence of the SSERKN methods compared with some well-known methods in the literature.     

Implicit space–time residual distribution method for unsteady laminar viscous flow

The class of multidimensional upwind residual distribution (RD) schemes has been developed in the past decades as an attractive alternative to the finite volume (FV) and finite element (FE) approaches. Although they have shown superior performances in the simulation of steady two-dimensional and three-dimensional inviscid and viscous flows, their extension to the simulation of unsteady flow fields is still a topic of intense research [ICCFD2, International Conference on Computational Fluid Dynamics 2, Sydney, Australia, 15–19 July 2002; M. Mezine, R. Abgrall, Upwind multidimensional residual schemes for steady and unsteady flows].Recently the space–time RD approach has been developed by several researchers [Int. J. Numer. Methods Fluids 40 (2002) 573; J. Comput. Phys. 188 (2003) 16; Á.G. Csı́k, Upwind residual distribution schemes for general hyperbolic conservation laws and application to ideal magnetohydrodynamics, PhD thesis, Katholieke Universiteit Leuven, 2002; J. Comput. Phys. 188 (2003) 16; R. Abgrall; M. Mezine, Construction of second order accurate monotone and stable residual distribution schemes for unsteady flow problems] which allows to perform second order accurate unsteady inviscid computations. In this paper we follow the work done in [Int. J. Numer. Methods Fluids 40 (2002) 573; Á.G. Csı́k, Upwind residual distribution schemes for general hyperbolic conservation laws and application to ideal magnetohydrodynamics, PhD thesis, Katholieke Universiteit Leuven, 2002]. In this approach the space–time domain is discretized and solved as a (d+1)-dimensional problem, where d is the number of space dimensions. In [Int. J. Numer. Methods Fluids 40 (2002) 573; Á.G. Csı́k, Upwind residual distribution schemes for general hyperbolic conservation laws and application to ideal magnetohydrodynamics, PhD thesis, Katholieke Universiteit Leuven, 2002] it is shown that thanks to the multidimensional upwinding of the RD method, the solution of the unsteady problem can be decoupled into sub-problems on space–time slabs composed of simplicial elements, allowing to obtain a true time marching procedure. Moreover, the method is implicit and unconditionally stable for arbitrary large time-steps if positive RD schemes are employed.We present further development of the space–time approach of [Int. J. Numer. Methods Fluids 40 (2002) 573; Á.G. Csı́k, Upwind residual distribution schemes for general hyperbolic conservation laws and application to ideal magnetohydrodynamics, PhD thesis, Katholieke Universiteit Leuven, 2002] by extending it to laminar viscous flow computations. A Petrov–Galerkin treatment of the viscous terms [Project Report 2002-06, von Karman Institute for Fluid Dynamics, Belgium, 2002; J. Dobeš, Implicit space–time method for laminar viscous flow], consistent with the space–time formulation has been investigated, implemented and tested. Second order accuracy in both space and time was observed on unstructured triangulation of the spatial domain.The solution is obtained at each time-step by solving an implicit non-linear system of equations. Here, following [Int. J. Numer. Methods Fluids 40 (2002) 573; Á.G. Csı́k, Upwind residual distribution schemes for general hyperbolic conservation laws and application to ideal magnetohydrodynamics, PhD thesis, Katholieke Universiteit Leuven, 2002], we formulate the solution of this system as a steady state problem in a pseudo-time variable. We discuss the efficiency of an explicit Euler forward pseudo-time integrator compared to the implicit Euler. When applied to viscous computation, the implicit method has shown speed-ups of more than a factor 50 in terms of computational time.     

Fast algorithms for high-order numerical methods for space-fractional diffusion equations

In this paper, fast numerical methods for solving space-fractional diffusion equations are studied in two stages. Firstly, a fast direct solver for an implicit finite difference scheme proposed by Hao et al. [A fourth-order approximation of fractional derivatives with its applications, J. Comput. Phys. 281 (2015), pp. 787–805], which is fourth-order accurate in space and second-order accurate in time, is developed based on a circulant-and-skew-circulant (CS) representation of Toeplitz matrix inversion. Secondly, boundary value method with spatial discretization of Hao et al. [A fourth-order approximation of fractional derivatives with its applications, J. Comput. Phys. 281 (2015), pp. 787–805] is adopted to produce a numerical solution with higher order accuracy in time. Particularly, a method with fourth-order accuracy in both space and time can be achieved. GMRES method is employed for solving the discretized linear system with two preconditioners. Based on the CS representation of Toeplitz matrix inversion, the two preconditioners can be applied efficiently, and the convergence rate of the preconditioned GMRES method is proven to be fast. Numerical examples are given to support the theoretical analysis.     

Space-time discretizing optimal DRP schemes for flow and wave propagation problems

The present paper reports constrained optimization of explicit Runge–Kutta (RK) schemes, coupled with optimal upwind compact scheme to achieve dispersion relation preservation (DRP) property for high performance computing. Essential ideas of optimization employed in arriving at the proposed time integration scheme are extension of the earlier work reported in Rajpoot et al. (J Comput Phys 2010;229:3623–51). This is in turn an application of the correct error evolution equation in Sengupta et al. (J Comput Phys 2007;226:1211–8). Resultant DRP scheme demonstrated the idea for explicit spatial central difference schemes. Present work is similar, extending it for near-spectral accuracy compact schemes. Practical utility of the developed method is demonstrated by solution of model problems and for flow problems by solving Navier–Stokes equation, some of which cannot be solved by conventional schemes, as the problem of rotary oscillation of cylinder.Developed method is calibrated with: (i) flow past a circular cylinder performing rotary oscillation at Re = 150 and (ii) flow inside a 2D lid-driven cavity (LDC) at Reynolds numbers of Re = 1000 and Re = 10,000. Quantitative and qualitative comparisons show excellent match for rotary oscillation cylinder cases with the experimental results of Thiria et al. (J Fluid Mech 2006;560:123–47). Results for LDC for Re = 1000 are compared with that in Botella & Peyret (Comp Fluids 1998;27:421–33) and results for Re = 10,000 are compared with recent published ones showing triangular vortex in the core.     

ADER schemes on unstructured meshes for nonconservative hyperbolic systems: Applications to geophysical flows

We develop a new family of well-balanced path-conservative quadrature-free one-step ADER finite volume and discontinuous Galerkin finite element schemes on unstructured meshes for the solution of hyperbolic partial differential equations with non-conservative products and stiff source terms. The fully discrete formulation is derived using the recently developed framework of explicit one-step P_NP_M schemes of arbitrary high order of accuracy in space and time for conservative hyperbolic systems [Dumbser M, Balsara D, Toro EF, Munz CD. A unified framework for the construction of one-step finite-volume and discontinuous Galerkin schemes. J Comput Phys 2008;227:8209–53]. The two key ingredients of our high order approach are: first, the high order accurate P_NP_M reconstruction operator on unstructured meshes, using the WENO strategy presented in [Dumbser M, Käser M, Titarev VA Toro EF. Quadrature-free non-oscillatory finite volume schemes on unstructured meshes for nonlinear hyperbolic systems. J Comput Phys 2007;226:204–43] to ensure monotonicity at discontinuities, and second, a local space–time Galerkin scheme to predict the evolution of the reconstructed polynomial data inside each element during one time step to obtain a high order accurate one-step time discretization. This approach is also able to deal with stiff source terms as shown in [Dumbser M, Enaux C, Toro EF. Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws. J Comput Phys 2008;227:3971–4001]. These two key ingredients are combined with the recently developed path-conservative methods of Parés [Parés C. Numerical methods for nonconservative hyperbolic systems: a theoretical framework. SIAM J Numer Anal 2006;44:300–21] and Castro et al. [Castro MJ, Gallardo JM, Parés C. High-order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. Applications to shallow-water systems. Math Comput 2006;75:1103–34] to treat the non-conservative products properly. We show applications of our method to the two-layer shallow water equations as well as applications to the recently published depth-averaged two-fluid flow model of Pitman and Le [Pitman EB, Le L. A two-fluid model for avalanche and debris flows. Philos Trans Roy Soc A 2005;363:1573–601].     

Characteristic-Based Schemes for Multi-Class Lighthill-Whitham-Richards Traffic Models

In this paper we provide the full spectral decomposition of the Multi-Class Lighthill Whitham Richards (MCLWR) traffic models described in (Wong et al. in Transp. Res. Part A 36:827–841, 2002; Benzoni-Gavage and Colombo in Eur. J. Appl. Math. 14:587–612, 2003). Even though the eigenvalues of these models can only be found numerically, the knowledge of the spectral structure allows the use of characteristic-based High Resolution Shock Capturing (HRSC) schemes. We compare the characteristic-based approach to the component-wise schemes used in (Zhang et al. in J. Comput. Phys. 191:639–659, 2003), and propose two strategies to minimize the oscillatory behavior that can be observed when using the component-wise approach.     

A Critical Study of the Compressible Lattice Boltzmann Methods for Riemann Problem

In this paper, the discrete velocity model proposed by Kataoka and Tsutahara (Phys. Rev. E 69(5):056702, 2004) for simulating inviscid flows is employed. Three approaches for improving the stability and the accuracy of this model, especially for high Mach numbers, are suggested and implemented in this research. First, the TVD scheme (Harten in J. Comput. Phys. 49:357?C393, 1983) is used for space discretization of the convective term in the Lattice Boltzmann equation. Next, the modified Lax-Wendroff with artificial viscosity is employed to increase the robustness of the method in supersonic flows. Finally, a combination of TVD and the 2nd order derivative of the distribution function is employed using a differentiable switch. It is found that the recent technique is a more suitable approach for a wide range of Mach numbers. Moreover, the WENO scheme for space discretization has been applied and compared with these newly applied methods.     

Vertex-Centered Linearity-Preserving Schemes for Nonlinear Parabolic Problems on Polygonal Grids

On arbitrary polygonal grids, a family of vertex-centered finite volume schemes are suggested for the numerical solution of the strongly nonlinear parabolic equations arising in radiation hydrodynamics and magnetohydrodynamics. We define the primary unknowns at the cell vertices and derive the schemes along the linearity-preserving approach. Since we adopt the same cell-centered diffusion coefficients as those in most existing finite volume schemes, it is required to introduce some auxiliary unknowns at the cell centers in the case of nonlinear diffusion coefficients. A second-order positivity-preserving algorithm is then suggested to interpolate these auxiliary unknowns via the primary ones. All the schemes lead to symmetric and positive definite linear systems and their stability can be rigorously analyzed under some standard and weak geometry assumptions. More interesting is that these vertex-centered schemes do not have the so-called numerical heat-barrier issue suffered by many existing cell-centered or hybrid schemes (Lipnikov et al. in J Comput Phys 305:111–126, 2016). Numerical experiments are also presented to show the efficiency and robustness of the schemes in simulating nonlinear parabolic problems.     

Guidelines for Poisson Solvers on Irregular Domains with Dirichlet Boundary Conditions Using the Ghost Fluid Method

We consider the variable coefficient Poisson equation with Dirichlet boundary conditions on irregular domains. We present numerical evidence for the accuracy of the solution and its gradients for different treatments at the interface using the Ghost Fluid Method for Poisson problems of Gibou et al. (J. Comput. Phys. 176:205–227, 2002; 202:577–601, 2005). This paper is therefore intended as a guide for those interested in using the GFM for Poisson-type problems (and by consequence diffusion-like problems and Stefan-type problems) by providing the pros and cons of the different choices for defining the ghost values and locating the interface. We found that in order to obtain second-order-accurate gradients, both a quadratic (or higher order) extrapolation for defining the ghost values and a quadratic (or higher order) interpolation for finding the interface location are required. In the case where the ghost values are defined by a linear extrapolation, the gradients of the solution converge slowly (at most first order in average) and the convergence rate oscillates, even when the interface location is defined by a quadratic interpolation. The same conclusions hold true for the combination of a quadratic extrapolation for the ghost cells and a linear interpolation. The solution is second-order accurate in all cases. Defining the ghost values with quadratic extrapolations leads to a non-symmetric linear system with a worse conditioning than that of the linear extrapolation case, for which the linear system is symmetric and better conditioned. We conclude that for problems where only the solution matters, the method described by Gibou, F., Fedkiw, R., Cheng, L.-T. and Kang, M. in (J. Comput. Phys. 176:205–227, 2002) is advantageous since the linear system that needs to be inverted is symmetric. In problems where the solution gradient is needed, such as in Stefan-type problems, higher order extrapolation schemes as described by Gibou, F. and Fedkiw, R. in (J. Comput. Phys. 202:577–601, 2005) are desirable.