Three-dimensional matched interface and boundary (MIB) method for treating geometric singularities
Introduction
The numerical solution of elliptic equations with discontinuous coefficients and singular sources has attracted much attention in the applied mathematics and scientific computing [5], [6], [7], [13], [16], [21], [26], [27], [28], [29], [31], [32], [34], [36], [41], [42], [43], [46], [61], [62], [63], [65], [66], [68] since Peskin’s pioneer work in 1977 [57]. This class of problems have demonstrated a growing impact to fluid dynamics [14], [18], [19], [30], [53], electromagnetics [23], [24], [33], material science [25], [27] and biologic systems [20], [44], [70], [75]. Due to the lack of regularity in the solution, ordinary numerical schemes, such as finite difference methods, finite element methods and spectral methods, cannot maintain the designed order of convergence and even diverge in the worst scenario. It is very difficult to construct highly accurate and efficient numerical methods for this class of problems. Apart from Peskin’s immersed boundary method (IBM) [22], [35], [57], [58], [59], a number of other elegant methods have been proposed in the literature. Among them, the ghost fluid method (GFM) proposed by Fedkiw, Osher and coworkers [15], [45] is a relatively simple and easy to use sharp interface approach. Finite element formulations [2], [8], [40], [44] are more suitable for weak solutions. An upwinding embedded boundary method was proposed by Cai and Deng [9] for dielectric media. A second-order finite volume based method was proposed Oevermann and Klein [54], which uses bilinear ansatz functions on Cartesian grids and compact stencils. A relevant, while quite distinct approach is the integral equation method for embedding complex geometry [49], [50], [51]. A major advance in the field was due to LeVeque and Li [37], who proposed a remarkable second-order sharp interface scheme, the immersed interface method (IIM). Their method has been improved for computational efficiency and matrix optimization over the past decade [1], [12], [38], [39], [62] and widely applied in practical problems. A new trend in this field is the construction of higher-order (interface) methods [21], [50], [63], [71], [73], [74] that are particularly desirable for problems involving both material interfaces and high frequency waves where conventional local adaptive refinement approaches do not work well. Typical examples are the interaction of turbulence and shock, and high frequency wave propagation in inhomogeneous media. Simple Cartesian grids are preferred in these situations because of the bypass of the mesh generation, better temporal stability and the availability of fast algebraic solvers. Recently, we have proposed a systematic higher-order method, the matched interface and boundary (MIB) method, for treating electromagnetic wave propagation and scattering in dielectric media [71]. In particular, we proposed a way to avoid the use of high order jump conditions, which require large stencils and are numerically unstable in higher-order schemes as analyzed in our work. We found that by repeatedly using only the lowest order jump conditions, an arbitrarily higher-order interface scheme can be constructed in principle [71]. This idea was applied to the fourth-order beam equation with free boundary conditions [72]. More recently, we have generalized the earlier MIB method for solving elliptic equations with curved interfaces [73], [74]. The MIB method makes use of simple Cartesian grids, standard finite difference schemes, lowest order physical jump conditions and fictitious domains. The physical jump conditions are enforced at each intersecting point of the interface and the mesh lines. Because of the use of the subgrid information at the intersecting points, the MIB method overcomes the staircase phenomena of the finite difference approximation for complex geometry and can be of arbitrarily higher-order convergence in principle when the lowest order physical jump conditions are repeatedly enforced. MIB schemes of up to 16th-order accuracy have been constructed for straight interfaces [71], [73]. For curved interfaces, fourth- and sixth-order MIB schemes have been demonstrated [73], [74]. Most recently, we have proposed an interpolation formulation of the MIB method that does not need the use of the fictitious values [74]. We have shown that our interpolation formulation of the MIB method is equivalent to our earlier fictitious domain formulation [74]. We have constructed a three-dimensional (3D) MIB scheme that is of second-order convergence for complex interfaces free of geometric singularities [75].
While much attention was focused on elliptic equations with smooth interfaces, little work has been done on the mathematical treatment of geometric singularities, such as sharp-edged, sharp-wedged and sharp-tipped interfaces. This class of problems are ubiquitous in science and engineering, such as wave-guides analysis [52], electromagnetic wave scattering and propagation [10], [55], [56], friction modeling [64], plasma-surface interaction [47] and turbulent-flow [4]. The presence of geometric singularities, such as cusps and self-intersecting surfaces, is a major obstacle to the accuracy, convergence and stability of the numerical solution of the Poisson–Boltzmann equation for the electrostatic analysis of biomolecules [70]. Technically, it is much more challenging to construct efficient numerical methods for this class of problems due to the lack of regularity in the interface. Because the gradient near the tip of a sharp-edged interface is not well defined in the classical sense, the aforementioned interface methods might not maintain the designed accuracy. A standard approach to these problems is the use of local adaptive approaches [52], [56], which exhibit reduced convergence rates and require dramatic mesh refinements in the vicinity of geometric singularities [11]. These approaches fail to work if the solution is highly oscillatory due to the so called pollution effect [3]. A remarkable result, due to Hou and Liu [26], was of 0.8th-order convergence obtained with a 2D finite element formulation for elliptic interface problems with geometric singularities. Their method is convenient for geometric singularity induced weak solutions. Recently, we have developed a 2D second-order MIB treatment of sharp-edged interfaces [69]. Unlike the Galerkin formulation, our MIB scheme does not specifically handle the geometric singularity induced weak solutions. A standard approach that multiplies the solution with an appropriate polynomial [56] was incorporated in our MIB method to deal with geometry induced weak solutions [69]. The MIB is found to be insensitive to the location of tips. The domain extension on the tip can be computed from any direction by extrapolation without affecting the result [69]. An explanation is that although gradients on the tip extrapolated from different directions contradict each other, different sets of interface conditions from different directions near the tip are consistent with each other for a given domain extension.
Another challenge that is encountered in interface schemes, particularly in higher-order ones, is the asymmetric discretization matrices which render convergence problems in solving the resulting linear algebraic equations. This difficulty is deteriorated by complex interface geometry—the matrices usually become more asymmetric under complex geometric constraints. Moreover, the higher convergence order, the larger matrix bandwidth will be for an interface scheme. Furthermore, discretization matrices are usually not diagonally dominant due to the interface treatment. This causes slow convergence as well. A discrete maximum preserving scheme was proposed to accelerate the convergence of the IIM [12]. In our earlier MIB method, a large number of iterations is required to solve the linear algebraic equations due to the asymmetric matrix.
The objective of the present work is threefold. First, we extend our earlier 2D higher-order MIB schemes [71], [73] to 3D. We report some first known sixth-order 3D results for elliptic problems with smooth interfaces. Such schemes are particularly efficient for problems involving both material interfaces and high frequency waves. Second, we generalize our earlier 2D second-order MIB method for geometric singularities [69] to 3D. Our earlier 3D MIB scheme does not maintain the designed second accuracy for molecular surfaces of proteins due to the presence of non-smooth interfaces [75]. In this work, we demonstrate the second-order convergence for arbitrarily complex interfaces with arbitrary geometric singularities. In fact, there is a critical acute angle restriction in our earlier 2D MIB scheme for sharp-edged interfaces [69]. A new algorithm, which makes use of two interface intersecting points for sharp-edged interfaces, is proposed to reserve this problem. We demonstrate that 3D fourth-order MIB schemes can be constructed for complex interfaces and interfaces with moderate geometric singularities. Finally, due to the constraint of complex interface geometry, it is impossible to create a high order MIB scheme with a symmetric discretization matrix. In this work, we have addressed this problem by making the MIB matrix optimally symmetric and banded. Here, by ‘optimally symmetric and banded’ we mean as symmetric and banded as possible under a given geometric constraint. The matrix optimization is achieved by an appropriate selection of auxiliary grid points. This optimization procedure dramatically improves the convergence property of the present MIB method. An essential feature of the proposed MIB method is that it locally simplifies a 2D or a 3D interface problem into 1D-like ones. As a result, the MIB method is able to systematically achieve higher-order convergence and to efficiently deal with geometric singularities. The proposed MIB algorithm, particularly its matrix acceleration and treatment of singularities, provides the crucial technical foundation to a new generation of Poisson–Boltzmann equation solvers for the electrostatic analysis of biomolecules [70].
In the next two sections, the theoretical formulation and the computational algorithm are given to the 3D MIB method for elliptic interface problems. Detailed consideration of topological variations is given, particularly for 3D geometric singularities. New algorithm is proposed for the sharp-edged interfaces and tips. The MIB matrix is made optimally symmetric and banded. A systematic procedure is introduced to construct 3D higher-order MIB schemes. The proposed MIB method is extensively validated by benchmark tests, such as the interfaces of missile, protein, flower and intersecting spheres. This paper ends with a conclusion to summarize the main points.
Section snippets
General setting of the MIB method for elliptic interface problems
Consider a given interface Γ which divides an open bounded domain into disjoint open subdomains, domain Ω− and domain Ω+, i.e., Ω = Ω− ∪ Ω+ ∪ Γ. We solve the following 3D elliptic equationwith variable coefficients β(x, y, z) and κ(x, y, z) which may admit jumps at the interface Γ. To make the problem well-posed, this equation should be solved with appropriate boundary conditions and two interface jump conditions
Results and discussion
In this section, we examine the convergence order, test the accuracy, validate the matrix optimization and demonstrate the capability of the proposed MIB II schemes. The impact of the proposed algorithm acceleration is studied with the molecular surfaces of a diatom and a protein. Molecular surfaces are generated with the MSMS program [60] with probe radius 1.4 and density 10. Reported CPU time is recorded on an AMD Turion 64™ MT-30 laptop with 1.6 GHz clock speed. The ability of the proposed
Conclusion
This paper introduces the three-dimensional (3D) higher-order matched interface and boundary (MIB) method for solving elliptic equations with discontinuous coefficients and geometric singularities. Our previous higher-order MIB methods in 1D [71] and 2D [73], [74] have been generalized to 3D. The present MIB method is of fourth-order convergence for complex interfaces with moderate geometric singularities, and of six-order convergence for smooth interfaces. Moreover, a recently proposed 2D
Acknowledgements
This work was supported in part by NSF Grants DMS-0616704 and IIS-0430987. The authors thank Shan Zhao, Yongcheng Zhou and Weihua Geng for useful discussions.
References (75)
- et al.
Numerical-simulation of turbulent-flow over surface-mounted obstacles with sharp edges and corners
J. Wind Engng. Indust. Aerodyn.
(1990) A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions
J. Comput. Phys.
(2004)- et al.
A fast solver for the Stokes equations with distributed forces in complex geometries
J. Comput. Phys.
(2004) - et al.
An upwinding embedded boundary method for Maxwell’s equations in media with material interfaces: 2D case
J. Comput. Phys.
(2003) - et al.
Three-dimensional elliptic solvers for interface problems and applications
J. Comput. Phys.
(2003) - et al.
A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)
J. Comput. Phys.
(1999) - et al.
A forth order accurate discretization for the Laplace and heat equations on arbitrary domains, with applications to the Stefan problem
J. Comput. Phys.
(2005) - et al.
On the order of accuracy of the immersed boundary method: higher order convergence rates for sufficiently smooth problems
J. Comput. Phys.
(2005) High-order accurate methods in time-domain computational electromagnetics. A review
Adv. Imaging Electron Phys.
(2003)- et al.
A numerical method for solving variable coefficient elliptic equation with interfaces
J. Comput. Phys.
(2005)