MIB method for elliptic equations with multi-material interfaces

Elliptic partial differential equations (PDEs) are widely used to model real-world problems. Due to the heterogeneous characteristics of many naturally occurring materials and man-made structures, devices, and equipments, one frequently needs to solve elliptic PDEs with discontinuous coefficients and singular sources. The development of high-order elliptic interface schemes has been an active research field for decades. However, challenges remain in the construction of high-order schemes and particularly, for nonsmooth interfaces, i.e., interfaces with geometric singularities. The challenge of geometric singularities is amplified when they are originated from two or more material interfaces joining together or crossing each other. High-order methods for elliptic equations with multi-material interfaces have not been reported in the literature to our knowledge. The present work develops matched interface and boundary (MIB) method based schemes for solving two-dimensional (2D) elliptic PDEs with geometric singularities of multi-material interfaces. A number of new MIB schemes are constructed to account for all possible topological variations due to two-material interfaces. The geometric singularities of three-material interfaces are significantly more difficult to handle. Three new MIB schemes are designed to handle a variety of geometric situations and topological variations, although not all of them. The performance of the proposed new MIB schemes is validated by numerical experiments with a wide range of coefficient contrasts, geometric singularities, and solution types. Extensive numerical studies confirm the designed second order accuracy of the MIB method for multi-material interfaces, including a case where the derivative of the solution diverges.     

Three-dimensional matched interface and boundary (MIB) method for treating geometric singularities

This paper reports the three-dimensional (3D) generalization of our previous 2D higher-order matched interface and boundary (MIB) method for solving elliptic equations with discontinuous coefficients and non-smooth interfaces. New MIB algorithms that make use of two sets of interface jump conditions are proposed to remove the critical acute angle constraint of our earlier MIB scheme for treating interfaces with sharp geometric singularities, such as sharp edges, sharp wedges and sharp tips. The resulting 3D MIB schemes are of second-order accuracy for arbitrarily complex interfaces with sharp geometric singularities, of fourth-order accuracy for complex interfaces with moderate geometric singularities, and of sixth-order accuracy for curved smooth interfaces. A systematical procedure is introduced to make the MIB matrix optimally symmetric and banded by appropriately choosing auxiliary grid points. Consequently, the new MIB linear algebraic equations can be solved with fewer number of iterations. The proposed MIB method makes use of Cartesian grids, standard finite difference schemes, lowest order interface jump conditions and fictitious values. The interface jump conditions are enforced at each intersecting point of the interface and mesh lines to overcome the staircase phenomena in finite difference approximation. While a pair of fictitious values are determined along a mesh at a time, an iterative procedure is proposed to determine all the required fictitious values for higher-order schemes by repeatedly using the lowest order jump conditions. A variety of MIB techniques are developed to overcome geometric constraints. The essential strategy of the MIB method is to locally reduce a 2D or a 3D interface problem into 1D-like ones. The proposed MIB method is extensively validated in terms of the order of accuracy, the speed of convergence, the number of iterations and CPU time. Numerical experiments are carried out to complex interfaces, including the molecular surfaces of a protein, a missile interface, and van der Waals surfaces of intersecting spheres.     

Time-Domain Numerical Solutions of Maxwell Interface Problems with Discontinuous Electromagnetic Waves

This paper is devoted to time domain numerical solutions of two-dimensional (2D) material interface problems governed by the transverse magnetic (TM) and transverse electric (TE) Maxwell's equations with discontinuous electromagnetic solutions. Due to the discontinuity in wave solutions across the interface, the usual numerical methods will converge slowly or even fail to converge. This calls for the development of advanced interface treatments for popular Maxwell solvers. We will investigate such interface treatments by considering two typical Maxwell solvers – one based on collocation formulation and the other based on Galerkin formulation. To restore the accuracy reduction of the collocation finite-difference time-domain (FDTD) algorithm near an interface, the physical jump conditions relating discontinuous wave solutions on both sides of the interface must be rigorously enforced. For this purpose, a novel matched interface and boundary (MIB) scheme is proposed in this work, in which new jump conditions are derived so that the discontinuous and staggered features of electric and magnetic field components can be accommodated. The resulting MIB time-domain (MIBTD) scheme satisfies the jump conditions locally and suppresses the staircase approximation errors completely over the Yee lattices. In the discontinuous Galerkin time-domain (DGTD) algorithm – a popular Galerkin Maxwell solver, a proper numerical flux can be designed to accurately capture the jumps in the electromagnetic waves across the interface and automatically preserves the discontinuity in the explicit time integration. The DGTD solution to Maxwell interface problems is explored in this work, by considering a nodal based high order discontinuous Galerkin method. In benchmark TM and TE tests with analytical solutions, both MIBTD and DGTD schemes achieve the second order of accuracy in solving circular interfaces. In comparison, the numerical convergence of the MIBTD method is slightly more uniform, while the DGTD method is more flexible and robust.     

Second-order Poisson Nernst-Planck solver for ion channel transport

The Poisson Nernst-Planck (PNP) theory is a simplified continuum model for a wide variety of chemical, physical and biological applications. Its ability of providing quantitative explanation and increasingly qualitative predictions of experimental measurements has earned itself much recognition in the research community. Numerous computational algorithms have been constructed for the solution of the PNP equations. However, in the realistic ion-channel context, no second order convergent PNP algorithm has ever been reported in the literature, due to many numerical obstacles, including discontinuous coefficients, singular charges, geometric singularities, and nonlinear couplings. The present work introduces a number of numerical algorithms to overcome the abovementioned numerical challenges and constructs the first second-order convergent PNP solver in the ion-channel context. First, a Dirichlet to Neumann mapping (DNM) algorithm is designed to alleviate the charge singularity due to the protein structure. Additionally, the matched interface and boundary (MIB) method is reformulated for solving the PNP equations. The MIB method systematically enforces the interface jump conditions and achieves the second order accuracy in the presence of complex geometry and geometric singularities of molecular surfaces. Moreover, two iterative schemes are utilized to deal with the coupled nonlinear equations. Furthermore, extensive and rigorous numerical validations are carried out over a number of geometries, including a sphere, two proteins and an ion channel, to examine the numerical accuracy and convergence order of the present numerical algorithms. Finally, application is considered to a real transmembrane protein, the Gramicidin A channel protein. The performance of the proposed numerical techniques is tested against a number of factors, including mesh sizes, diffusion coefficient profiles, iterative schemes, ion concentrations, and applied voltages. Numerical predictions are compared with experimental measurements.     

The adaptive immersed interface finite element method for elliptic and Maxwell interface problems

We propose self-adaptive finite element methods with error control for solving elliptic and electromagnetic problems with discontinuous coefficients. The meshes in the methods do not need to fit the interfaces. New error indicators are introduced to control the error due to non-body-fitted meshes. Flexible h-adaptive strategies are developed, which can be systematically extended to a large class of interface problems. Extensive numerical experiments are performed to support the theoretical results and to show the competitive behavior of the adaptive algorithm even for interfaces involving corner or tip singularities.     

Modelling atmospheric flows with adaptive moving meshes

An anelastic atmospheric flow solver has been developed that combines semi-implicit non-oscillatory forward-in-time numerics with a solution-adaptive mesh capability. A key feature of the solver is the unification of a mesh adaptation apparatus, based on moving mesh partial differential equations (PDEs), with the rigorous formulation of the governing anelastic PDEs in generalised time-dependent curvilinear coordinates. The solver development includes an enhancement of the flux-form multidimensional positive definite advection transport algorithm (MPDATA) — employed in the integration of the underlying anelastic PDEs — that ensures full compatibility with mass continuity under moving meshes. In addition, to satisfy the geometric conservation law (GCL) tensor identity under general moving meshes, a diagnostic approach is proposed based on the treatment of the GCL as an elliptic problem. The benefits of the solution-adaptive moving mesh technique for the simulation of multiscale atmospheric flows are demonstrated. The developed solver is verified for two idealised flow problems with distinct levels of complexity: passive scalar advection in a prescribed deformational flow, and the life cycle of a large-scale atmospheric baroclinic wave instability showing fine-scale phenomena of fronts and internal gravity waves.     

Three-dimensional elliptic solvers for interface problems and applications

Second-order accurate elliptic solvers using Cartesian grids are presented for three-dimensional interface problems in which the coefficients, the source term, the solution and its normal flux may be discontinuous across an interface. One of our methods is designed for general interface problems with variable but discontinuous coefficient. The scheme preserves the discrete maximum principle using constrained optimization techniques. An algebraic multigrid solver is applied to solve the discrete system. The second method is designed for interface problems with piecewise constant coefficient. The method is based on the fast immersed interface method and a fast 3D Poisson solver. The second method has been modified to solve Helmholtz/Poisson equations on irregular domains. An application of our method to an inverse interface problem of shape identification is also presented. In this application, the level set method is applied to find the unknown surface iteratively.     

A Cartesian grid embedded boundary method for solving the Poisson and heat equations with discontinuous coefficients in three dimensions

We present a method for solving Poisson and heat equations with discontinuous coefficients in two- and three-dimensions. It uses a Cartesian cut-cell/embedded boundary method to represent the interface between materials, as described in Johansen and Colella (1998). Matching conditions across the interface are enforced using an approximation to fluxes at the boundary. Overall second order accuracy is achieved, as indicated by an array of tests using non-trivial interface geometries. Both the elliptic and heat solvers are shown to remain stable and efficient for material coefficient contrasts up to 10⁶, thanks in part to the use of geometric multigrid. A test of accuracy when adaptive mesh refinement capabilities are utilized is also performed. An example problem relevant to nuclear reactor core simulation is presented, demonstrating the ability of the method to solve problems with realistic physical parameters.     

A domain decomposition method for solving the three-dimensional time-harmonic Maxwell equations discretized by discontinuous Galerkin methods

We present here a domain decomposition method for solving the three-dimensional time-harmonic Maxwell equations discretized by a discontinuous Galerkin method. In order to allow the treatment of irregularly shaped geometries, the discontinuous Galerkin method is formulated on unstructured tetrahedral meshes. The domain decomposition strategy takes the form of a Schwarz-type algorithm where a continuity condition on the incoming characteristic variables is imposed at the interfaces between neighboring subdomains. A multifrontal sparse direct solver is used at the subdomain level. The resulting domain decomposition strategy can be viewed as a hybrid iterative/direct solution method for the large, sparse and complex coefficients algebraic system resulting from the discretization of the time-harmonic Maxwell equations by a discontinuous Galerkin method.     

On the elliptic mesh generation in domains containing multiple inclusions and undergoing large deformations

We present an improved method to generate a sequence of structured meshes even when the physical domain contains deforming inclusions. This method belongs to the class of Arbitrary Lagrangian–Eulerian (ALE) methods for solving moving boundary problems. Its tools are either (a) separate mappings of the domain boundaries and enforcing the node distribution on lines emanating from singular points or (b) domain decomposition and separate mappings of each subdomain using suitable coordinate systems. The latter is shown to be more versatile and general. In both cases a set of elliptic equations is used to generate the grid extending in this way the method advanced by Dimakopoulos and Tsamopoulos [Y. Dimakopoulos, J.A. Tsamopoulos, A quasi-elliptic transformation for moving boundary problems with large anisotropic deformations, J. Comput. Phys. 192 (2003) 494–522]. We shall present examples where this earlier method and all other mesh generating methods which are based on a conformal mapping or solving a quasi-elliptic set of PDEs fail to produce an acceptable mesh and accurate solutions in such geometries. Furthermore, in contrast to other methods, appropriate boundary conditions and constraints such as, orthogonality of specific mesh lines and prespecified node distributions on them, can be easily implemented along a specific part of the domain or its boundary. Hence, no attractive terms at specific corners or singular points are needed. To increase the mesh resolution around the moving interfaces while keeping low the memory requirements and the computational time, a local mesh refinement technique has been incorporated as well. The method is demonstrated in two challenging examples where no remeshing is required in spite of the large domain deformations. In the first one, the transient growth of two bubbles embedded in a viscoelastic filament undergoing stretching in the axial direction is examined, while in the second one the linear and non-linear dynamics of two bubbles in a viscous medium are determined in an acoustic field. The large elasticity of the filament in the first case or the large inertia in the second case coupled with the externally induced large deformations of the liquid domain requires the accurate calculation which is achieved by the method we propose herein. The governing equations are solved using the finite element/Galerkin method with appropriate modifications to solve the hyperbolic constitutive equation of a viscoelastic fluid. These are coupled with an implicit Euler method for time integration or with Arnoldi’s algorithm for normal mode analysis.     

流体力学方程的间断有限元方法

在二维区域三角形网格上应用一阶、二阶和三阶精度间断有限元方法,对流体力学方程和方程组进行了数值模拟.计算结果与差分方法计算结果比较,认为间断有限元方法在求解复杂边界条件和区域问题上有一定的优势.     

An Error Indicator Monitor Function for an r-Adaptive Finite-Element Method

An r-adaptive finite-element method based on moving-mesh partial differential equations (PDEs) and an error indicator is presented. The error indicator is obtained by applying a technique developed by Bank and Weiser to elliptic equations which result in this case from temporal discretization of the underlying physical PDEs on moving meshes. The construction of the monitor function based on the error indicator is discussed. Numerical results obtained with the current method and the commonly used method based on solution gradients are presented and analyzed for several examples.     

Locally implicit discontinuous Galerkin method for time domain electromagnetics

In the recent years, there has been an increasing interest in discontinuous Galerkin time domain (DGTD) methods for the solution of the unsteady Maxwell equations modeling electromagnetic wave propagation. One of the main features of DGTD methods is their ability to deal with unstructured meshes which are particularly well suited to the discretization of the geometrical details and heterogeneous media that characterize realistic propagation problems. Such DGTD methods most often rely on explicit time integration schemes and lead to block diagonal mass matrices. However, explicit DGTD methods are also constrained by a stability condition that can be very restrictive on highly refined meshes and when the local approximation relies on high order polynomial interpolation. An implicit time integration scheme is a natural way to obtain a time domain method which is unconditionally stable but at the expense of the inversion of a global linear system at each time step. A more viable approach consists of applying an implicit time integration scheme locally in the refined regions of the mesh while preserving an explicit time scheme in the complementary part, resulting in an hybrid explicit–implicit (or locally implicit) time integration strategy. In this paper, we report on our recent efforts towards the development of such a hybrid explicit–implicit DGTD method for solving the time domain Maxwell equations on unstructured simplicial meshes. Numerical experiments for 3D propagation problems in homogeneous and heterogeneous media illustrate the possibilities of the method for simulations involving locally refined meshes.     

Covolume-upwind finite volume approximations for linear elliptic partial differential equations

In this paper, we study covolume-upwind finite volume methods on rectangular meshes for solving linear elliptic partial differential equations with mixed boundary conditions. To avoid non-physical numerical oscillations for convection-dominated problems, nonstandard control volumes (covolumes) are generated based on local Peclet’s numbers and the upwind principle for finite volume approximations. Two types of discretization schemes with mass lumping are developed with use of bilinear or biquadratic basis functions as the trial space respectively. Some stability analyses of the schemes are presented for the model problem with constant coefficients. Various examples are also carried out to numerically demonstrate stability and optimal convergence of the proposed methods.     

A new anisotropic mesh adaptation method based upon hierarchical a posteriori error estimates

A new anisotropic mesh adaptation strategy for finite element solution of elliptic differential equations is presented. It generates anisotropic adaptive meshes as quasi-uniform ones in some metric space, with the metric tensor being computed based on hierarchical a posteriori error estimates. A global hierarchical error estimate is employed in this study to obtain reliable directional information of the solution. Instead of solving the global error problem exactly, which is costly in general, we solve it iteratively using the symmetric Gauß–Seidel method. Numerical results show that a few GS iterations are sufficient for obtaining a reasonably good approximation to the error for use in anisotropic mesh adaptation. The new method is compared with several strategies using local error estimators or recovered Hessians. Numerical results are presented for a selection of test examples and a mathematical model for heat conduction in a thermal battery with large orthotropic jumps in the material coefficients.     

High order matched interface and boundary methods for the Helmholtz equation in media with arbitrarily curved interfaces

This work overcomes the difficulty of the previous matched interface and boundary (MIB) method in dealing with interfaces with non-constant curvatures for optical waveguide analysis. This difficulty is essentially bypassed by avoiding the use of local cylindrical coordinates in the improved MIB method. Instead, novel jump conditions are derived along global Cartesian directions for the transverse magnetic field components. Effective interface treatments are proposed to rigorously impose jump conditions across arbitrarily curved interfaces based on a simple Cartesian grid. Even though each field component satisfies the scalar Helmholtz equation, the enforcement of jump conditions couples two transverse magnetic field components, so that the resulting MIB method is a full-vectorial approach for the modal analysis of optical waveguides. The numerical performance of the proposed MIB method is investigated by considering interface problems with both constant and general curvatures. The MIB method is shown to be able to deliver a fourth order of accuracy in all cases, even when a high frequency solution is involved.     

A spectral fictitious domain method with internal forcing for solving elliptic PDEs

A fictitious domain method is presented for solving elliptic partial differential equations using Galerkin spectral approximation. The fictitious domain approach consists in immersing the original domain into a larger and geometrically simpler one in order to avoid the use of boundary fitted or unstructured meshes. In the present study, boundary constraints are enforced using Lagrange multipliers and the novel aspect is that the Lagrange multipliers are associated with smooth forcing functions, compactly supported inside the fictitious domain. This allows the accuracy of the spectral method to be preserved, unlike the classical discrete Lagrange multipliers method, in which the forcing is defined on the boundaries. In order to have a robust and efficient method, equations for the Lagrange multipliers are solved directly with an influence matrix technique. Using a Fourier–Chebyshev approximation, the high-order accuracy of the method is demonstrated on one- and two-dimensional elliptic problems of second- and fourth-order. The principle of the method is general and can be applied to solve elliptic problems using any high order variational approximation.     

Reduced basis finite element heterogeneous multiscale method for high-order discretizations of elliptic homogenization problems

A new finite element method for the efficient discretization of elliptic homogenization problems is proposed. These problems, characterized by data varying over a wide range of scales cannot be easily solved by classical numerical methods that need mesh resolution down to the finest scales and multiscale methods capable of capturing the large scale components of the solution on macroscopic meshes are needed. Recently, the finite element heterogeneous multiscale method (FE-HMM) has been proposed for such problems, based on a macroscopic solver with effective data recovered from the solution of micro problems on sampling domains at quadrature points of a macroscopic mesh. Departing from the approach used in the FE-HMM, we show that interpolation techniques based on the reduced basis methodology (an offline-online strategy) allow one to design an efficient numerical method relying only on a small number of accurately computed micro solutions. This new method, called the reduced basis finite element heterogeneous multiscale method (RB-FE-HMM) is significantly more efficient than the FE-HMM for high order macroscopic discretizations and for three-dimensional problems, when the repeated computation of micro problems over the whole computational domain is expensive. A priori error estimates of the RB-FE-HMM are derived. Numerical computations for two and three dimensional problems illustrate the applicability and efficiency of the numerical method.     

Krylov implicit integration factor methods for spatial discretization on high dimensional unstructured meshes: Application to discontinuous Galerkin methods

Integration factor methods are a class of “exactly linear part” time discretization methods. In [Q. Nie, Y.-T. Zhang, R. Zhao, Efficient semi-implicit schemes for stiff systems, Journal of Computational Physics, 214 (2006) 521–537], a class of efficient implicit integration factor (IIF) methods were developed for solving systems with both stiff linear and nonlinear terms, arising from spatial discretization of time-dependent partial differential equations (PDEs) with linear high order terms and stiff lower order nonlinear terms. The tremendous challenge in applying IIF temporal discretization for PDEs on high spatial dimensions is how to evaluate the matrix exponential operator efficiently. For spatial discretization on unstructured meshes to solve PDEs on complex geometrical domains, how to efficiently apply the IIF temporal discretization was open. In this paper, we solve this problem by applying the Krylov subspace approximations to the matrix exponential operator. Then we apply this novel time discretization technique to discontinuous Galerkin (DG) methods on unstructured meshes for solving reaction–diffusion equations. Numerical examples are shown to demonstrate the accuracy, efficiency and robustness of the method in resolving the stiffness of the DG spatial operator for reaction–diffusion PDEs. Application of the method to a mathematical model in pattern formation during zebrafish embryo development shall be shown.     

Coefficient Jump-Independent Approximation of the Conforming and Nonconforming Finite Element Solutions

A counterexample is constructed. It confirms that the error of conforming finite element solution is proportional to the coefficient jump, when solving interface elliptic equations. The Scott-Zhang operator is applied to a nonconforming finite element. It is shown that the nonconforming finite element provides the optimal order approximation in interpolation, in $L^2$-projection, and in solving elliptic differential equation, independent of the coefficient jump in the elliptic differential equation. Numerical tests confirm the theoretical finding.