A multigrid based finite difference method for solving parabolic interface problem

1.   Introduction

The paper is focused on the following two dimensional (2D) parabolic interface problem

with Dirichlet boundary conditions imposed along the boundary $\partial \Omega$ of a given rectangular domain. The whole domain is separated by an embedded interface into two subdomains such that $\Gamma = \Omega^+ \cap \Omega^-,$ and $\Omega = \Omega^+ \cup \Omega^-.$ In our study, we specify $\Omega^+$ as outer domain, and $\Omega^-$ as inner domain. The diffusion coefficient $\beta$ is defined to be a piecewise constant in different subdomains. Across the interface $\Gamma$, the function values of $u$ from the two subdomains are associated by the following jump conditions

where superscripts $+$ and $-$ on $u$ and $\beta$ signify the corresponding subdomains, $\vec{n}$ indicates the normal direction pointing from $\Omega^-$ to $\Omega^+$, and the jump values $\phi$ and $\psi$ are allowed to be temporally and spatially dependent to account for the most general scenario.

The parabolic interface problem modeled by (1) - (3) plays important roles in multiple physical and engineering applications. The conductive heat transfer process over composite media is widely investigated through this model when the interface conditions are homogeneous, i.e., $[\![ u ]\!] = 0$, and $[\![ \beta u_n ]\!] = 0.$ Due to the discontinuity of diffusion coefficient across the interface, the analytical solution of the parabolic interface problem is not readily available even in the case of homogeneous interface jump conditions. Hence numerical approaches are in need for such problems. However, standard numerical procedures cannot deliver accurate solutions due to the assumption that solution is smooth over the whole domain.

Various specially designed numerical schemes with the jump conditions being incorporated into discretization have been introduced in the literature for solving parabolic interface problems. These include finite element methods and finite volume methods based on body-fitted interface treatments [3,8,32,33,34], finite difference methods based on Cartesian grids [2,5,6,17,18,19,20,22,24,30,36,40], and immersed finite element methods based on Cartesian meshes [25,39]. For instance, the immersed interface method (IIM) is one of the most successful finite difference methods in solving parabolic interface problems [5,6,18,19,20]. By avoiding mesh generation, generalized Taylor expansions are adopted to locally modify the finite difference weights so that the accuracy loss across the interface can be recovered.

Time integration is also crucial to the numerical solution of parabolic interface problems. Explicit time stepping would require a severe stability restriction with the time step size $\Delta t<Ch^2/\beta_{max}$, where $h$ is the spatial grid size and $\beta_{max} = \text{max}_{x\in \Omega}\beta(x)$. The implicit time integration is thus preferred in many applications involving long-time integration for the steady state solution. However, in implicit time stepping, a linear system, which involves all spatial degree of freedom $N$ and is obtained from discretizing a Helmholtz type equation, has to be solved in each time step. With a generic iterative solver, the computational complexity could be $O(N^2)$ in each time step. This leads to tremendous amount of computing resources, especially when the parabolic interface problem is generalized to three dimensions. This calls for the development of fast algorithms for solving such a linear system with the complexity being $O(N)$ or $O(N\log N)$. In scientific computing, some possible candidates with such efficiency could be the fast Fourier transform (FFT), multigrid, or alternating directional implicit (ADI) methods. Nevertheless, the FFT is not readily applicable to such parabolic interface problem due to the non-constant coefficient of the reaction term in the Helmholtz equation. Therefore, ADI and multigrid are two popular techniques for designing fast algorithms for parabolic partial differential equations (PDEs).

The core philosophy of ADI methods [12,11,29] is to reduce a multidimensional system from discretizing parabolic equation to a sequence of independent one-dimensional (1D) sub-systems of tridiagonal structure, which can be solved efficiently with the Thomas algorithm [31]. Compared with generic iterative solvers, ADI methods can be regarded as an exact solver due to the fact that computation will be completed within a fixed number of steps. Recently, ADI schemes have been developed to cope with parabolic interface problems with desired ADI efficiency maintained. Li and Mayo [22] introduced the first ADI algorithm with help of the IIM to correct local differences near interface when facing a special interface jump with $[\![\beta]\!] = 0$ across the interface. The further studies of such IIM scheme can be found in [26,27]. Zhao [40] proposed the first ADI algorithm to tackle the general jump conditions with $[\![\beta]\!]\ne 0$, by adopting matched interface and boundary (MIB) method [41,42] for handling interfaces. A tensor product approach is applied to decompose the jump conditions into Cartesian and tangential directions. The Cartesian jump quantities can be used to correct 1D finite difference discretization, while the tangential jump is approximated by using function values in the previous time step. A modified Thomas algorithm is introduced to accomplish the classic computational complexity $O(N)$ of the ADI scheme. The same MIB-based philosophy is further generalized to three dimensional parabolic interface problem with a non-orthogonal local coordinate for a simpler implementation [36]. These MIB based ADI schemes [40,36] achieve second order accuracy in space, but attain only a first order accuracy in time, due to tangential approximations. More recently, a multiscale ADI algorithm [20] has been proposed with aid of the IIM for general interface conditions. Two types of auxiliary variables on the interface are introduced for the augmented system in order to incorporate the ADI framework. Second order accuracy is fulfilled for both spatial and temporal discretization with a time step size constraint $\Delta t\leq Ch.$ Most recently, the ghost-fluid method has also been combined with the ADI to solve parabolic interface problem with symmetric and tridiagonal finite difference matrix obtained [23]. The metric of the algorithm lies in its superior stability despite first order accuracy in space.

The multigrid method is another popular fast solver used in solving PDEs with discontinuous coefficients. Based on structured or unstructured meshes, approaches for constructing interpolation, smoothing operators and coarse gird points selection would require special designs [7]. Algebraic multigrid (AMG) [28] selects coarse grid points in a pure algebraic sense to define interpolation without interface geometry and dimensionality needed in the construction. Yet it exploits the discontinuous coefficients and geometry in the matrix-dependent interpolation. A great amount of storage is usually needed for its coarse grid operator matrices generated by Galerkin process. Recently, the AMG scheme has been successfully applied in the immersed finite element method for solving both stationary and moving interface problems [13]. Geometric multigrid method is the other class of multigrid methods, which includes interface geometry in the process of discretization to the PDE with related boundary and interface jump conditions. In [2], Adams and Li proposed a 2D IIM multigrid method, which preserves maximum principle. Black box multigrid interpolation is adopted away from the interface, while interpolation weights are derived with help of Taylor expansion for the grid points adjacent to the interface. Later on, an improved IIM based multigrid method [1] has been developed with interpolation and restriction operators being modified such that coarse-grid matrices are M-matrices. The number of $V$-cycles does not change with mesh refinement, and increases slightly with the ratio of discontinuous coefficients growing. In [35], a multigrid method is proposed, which can maintain a uniform convergence with respect to mesh size and the high interface jump ratio. It can be extended to three dimensions, and can solve multiple interface problem efficiently. Besides, it can avoid the issue of large storage of coarse-grid matrix. Extending the multigrid approach for solving continuous coefficient boundary value problem in [9], Coco further designed a multigrid method for solving discontinuous coefficient interface problem [10], which can even be further generalized to matrix coefficient problems. Its multigrid efficiency is guaranteed with the convergence factor not depending on the coefficient ratio. In spite of the great success achieved by the aforementioned multigrid methods, a new multigrid method has to be designed again when a new interface algorithm is under consideration. Because each interface algorithm has a different procedure to capture boundary or interface conditions, sophisticated treatments are demanded to define interpolation with knowledge of geometry from the interface and to define restriction and coarse grid operators.

In this paper, we propose a new geometric multigrid method based on the augmented matched interface and boundary (AMIB) method for solving parabolic interface problems. The AMIB method was originally designed to solve elliptic interface problems [14,15,16]. With a FFT acceleration, the AMIB method can deliver fourth order accuracy in treating both interfaces [16] and boundaries [15]. Nevertheless, the AMIB method has never been coupled with a multigrid method before. In this study, due to the augmented approach adopted, the finite difference interface treatment and the multigrid solution procedure can be decoupled in the proposed AMIB method. This considerably simplifies the interpolation and restriction process for the coarse grid operators in the multigrid method. Hence standard multigrid components can be applied in such augmented framework with the advantages including that storage for coarse grid matrix is avoided and the implementation is fairly simple. Second order corrected central differences are applied for spatial accuracy to compensate accuracy loss across the interface, and the Crank-Nicolson scheme is used for temporal discretization. This leads to second order accuracy in both space and time, while the multigrid approach considerably accelerates the computational speed.

The rest of the paper is organized as follows. Section 2 focuses on discussion of the several key aspects of the proposed efficient algorithm. In section 3, numerical experiments are implemented to justify the effectiveness of our method in solving parabolic interface problems. Comparisons to other methods are carried out to show the efficiency of our method. A summary and acknowledgment are included in the end of the paper.

2.   Theory and algorithm

In this work, we focus on the piecewise constant coefficient parabolic interface problem. We first rewrite the original problem $(1)$ into the following form

Subject to the same initial, boundary, and interface conditions, the solution to (4) is equivalent to that of (1). We partition the given rectangular domain $\Omega = [a, b]\times[c, d]$ into to $n_x$ and $n_y$ equally spaced intervals in $x-$ and $y-$ directions respectively. The mesh size in each direction can then be defined as $h_x = (b-a)/n_x$ and $h_y = (d-c)/n_y$. For simplicity, we assume that $h_x = h_y = h.$ The grid nodes are hence defined as

Assume that the interface $\Gamma$ is represented by a level set function

with $\varphi(x, y) > 0$ in $\Omega^+$ and $\varphi(x, y) < 0$ in $\Omega^-$. Two functions are defined at each grid point

The subscript $(i, j)$ indicates the function value of $\varphi$ at $(x_i, y_j)$. It is similarly defined for other variables. If $\varphi_{i, j}^{min}\varphi_{i, j}^{max}>0$, then the grid point $(x_i, y_j)$ is called a regular point, otherwise an irregular point.

Moreover, a uniform time increment $\triangle t$ is adopted for temporal discretization. Let $t_n = n\cdot \triangle t$ represent the $n$th time step in the temporal discretization. For ease of exposition, we use notation $u_{i, j}^n$ to denote the function value on grid $(x_i, y_j)$ at time $t_n$.

2.1.   Temporal discretization

For temporal discretization, the Crank-Nicolson scheme is adopted

which has the equivalent form below after reorganizing the terms in (5)

The above Crank-Nicolson scheme admits a second order temporal truncation error. Note that the Laplacian terms at time step $t = t^{n+1}$ and $t = t^n$ need spatial discretization such that the time stepping process (6) can proceed for the discrete solutions at certain time instant. The spatial discretization details are fully demonstrated in following sections.

2.2.   Spatial discretization

Ignoring the time dependence, we concern ourselves with second order spatial discretization to the Laplacian operator

We may focus the discussion on $\frac{\partial^2u}{\partial x^2}$ due to the tensor-product decomposition to be employed for discretizing the Laplacian operator. The standard second order central difference can be utilized to approximate the partial derivative $u_{xx}$ with a truncation error $O(h^2)$:

The above formulation (7) is applicable to the case of regular points, while approximation in irregular points needs different treatment due to the fact that function loses regularity across the interface. Corrected differences are employed below [37]:

Jump-corrected difference. Suppose $u\in C^4[x_i-h, \alpha)\cap C^4(\alpha, x_{i+1}+h]$ and $x_i\in \Omega^{-}$ and $x_{i+1}\in \Omega^{+}$, where derivatives extend continuously up to $\alpha$. The following approximations hold to $O(h^2)$ when $K = 3$ :

where $h^{-} = x_i-\alpha$, $h^{+} = x_{i+1}-\alpha$ with $x_i\leq\alpha<x_{i+1}$, and $[u^{(k)}]_{\alpha} = \lim\limits_{x\rightarrow \alpha^{+}}u^{(k)}(x)-\lim\limits_{x\rightarrow \alpha^{-}}u^{(k)}(x).$

The global regularity for $u$ will be assumed to be piecewise $C^4$ continuous, as required locally for jump-corrected differences. Note that the leading error term in (7) involves the fourth order derivative $u_{xxxx}$. So $C^4$ continuity is needed in each subdomain for maintaining the second order accuracy everywhere, except near the interface. In corrected differences, the Cartesian derivative jumps of up to order $K = 2$ is sufficient to guarantee a global second order accuracy despite being first order locally [14].

The corrected differences in Eqns.(8) and (9) will reduce to Eq. (7) with the introduced jump terms vanishing if no interface is met. Such jump quantities are compensating for the function discontinuity across the interface. Since the derivative jumps in the correction terms are not analytically available, numerical approximations are needed for the corrected differences to be applied at irregular points. In following subsection, a systematical procedure is discussed for reconstruction of the derivative jumps with aid of fictitious values generated by the matched interface and boundary (MIB) method [42].

2.3.   Reconstructing the Cartesian derivative jumps

The derivative jumps in the above corrected differences are defined as

where $k = 0, 1, 2$ indicates the order of derivative. We will construct an effective numerical approximation to such derivative jumps except for the zeroth order one, which is known in the interface condition (2). It is natural to approximate the one-side derivative jumps in Eq.(10) using polynomials from each side. Since the derivative order in Eq.(10) is up to two, a polynomial of degree two on each side of the interface is sufficient to guarantee the local first order accuracy in the corrected difference with each order derivative obtained via differentiation on the polynomial.

For the sake of accuracy and stability, it is better to approximate each one-sided derivative limit with central or pseudo-central difference. As shown in Fig. 1, if one layer of fictitious values is available on each side of the interface $x = \alpha$, we can combine those with two real function values on the other side of the interface to yield the desired two Lagrange polynomials. The derivative jumps can be approximated as below:

where the linear combinations in the two parenthesis denote the derivatives of the two constructed Lagrange polynomials for approximating one-sided derivative limits at $x = \alpha$.

It is possible that only one real function value is available on one side of the interface due to sharp topology change of the interface. To deal with such a corner case [14], a fictitious value is employed at where the needed real value is missing. Just consider the derivative jumps at $x = \alpha_2$ in the two scenarios shown in Fig. 2. The fictitious value at $x = x_{i-1}$ is used for the case in the right chart in comparison with the case in the left chart where the real value at $x = x_{i-1}$ is adopted.

Likewise, the discretization for the $y-$direction derivative jump can be formulated when the interface intersects with the $x$ grid line. To apply the approximation (11), we need to construct fictitious values at all irregular points. In the following subsection, the MIB method is employed to generate the needed fictitious values.

2.4.   Fictitious values formulation

The jump conditions $(2)$ and $(3)$ are given on the interface $\Gamma$. One more interface condition can be analytically derived by differentiating Eq. (2) along the tangential direction $\tau$ of the interface as below,

Consider a point $(x^*, y^*)$ on the interface. Define $\theta$ as the angle between positive $x-$direction and normal direction. Then the tangential and normal direction could be expressed as $\vec{\tau} = (-\sin\theta, \cos\theta)$ and $\vec{n} = (\cos\theta, \sin\theta)$, respectively. See Fig. 3 in case that $(x^*, y^*)$ is located on a $y$ grid line. At $(x^*, y^*)$, the three known interface conditions can be reinterpreted in the Cartesian form [42],

In this way, jump conditions from tangential or normal direction are transformed into $x-$or $y-$direction involving four partial derivatives $u_x^+, u_x^-, u_y^+,$ and $u_y^-$. The following formulations (16) or (17) can be obtained by eliminating $u_y^-$ or $u_x^-$ from the four quantities.

If $u_y^-$ is eliminated, the following equations can be derived

where $C_x^{+} = \beta^{+}\cos \theta+\beta^{-}\tan \theta \sin \theta$, $C_x^{-} = \beta^{-}\cos \theta+\beta^{-}\tan\theta \sin\theta$ and $C_y^{+} = \beta^{+}\sin \theta-\beta^{-}\sin \theta.$

If $u_x^-$ is removed, we can obtain that

where $D_x^{+} = (\beta^{+}-\beta^{-})\cos\theta, D_y^{+} = \beta^{-}\cos\theta \cot\theta+\beta^{+}\sin\theta$ and $D_y^{-} = \beta^{-}(\cos\theta \cot\theta +\sin\theta)$.

Eqns. (16) and (17) allow us to reduce the 2D jump conditions to 1D ones along Cartesian direction provided that $u_y^{+}$ and $u_x^{+}$ can be determined, respectively. Suppose the interface intersects with $y = y_j$ at $(x_o, y_j)$ as shown in Fig. 3. To generate a pair of fictitious values $\hat{u}_{i, j}$ and $\hat{u}_{i+1, j}$ at points $P_1$ and $P_2$ on each side of the interface, discretization is carried out on Eqn.(16) as below,

with the following vector representations

Here $U^+, U^-$ indicate the involved function and fictitious values and the four vectors $W_0^+$, $W_0^-$, $W_1^+$, $W_1^-$ represent the Lagrange interpolation weights when discretizing (16) along $x-$direction. Note that fictitious values $\hat{u}_{i, j}$ and $\hat{u}_{i+1, j}$ in $U^+$ and $U^{-}$ participate in the process of approximation to the $x-$direction interpolation. Analogously, $U_{o}$ and $P^+$ are the needed auxiliary function values and corresponding weights to approximate $u_y^{+}$. The superscripts $-$ and $+$ in the above vector notations or the elements in each vector signify the $\Omega^{-}$ and $\Omega^{+}$ domain. Beside, the subscript $0$ or $1$ indicates the zeroth or first order derivative. The index $i$ or its variance tells the information located at $x = x_i$. Note that three auxiliary function values in $U_{o}$ are at points $(x_{o}, y_{j-1}), (x_{o}, y_j)$, $(x_{o}, y_{j+1})$, see Fig 3. Additionally, each of these three auxiliary values are interpolated or extrapolated by three function values in $\Omega^{+}$. For example, in Fig. 3, $u_{o, j+1}$ is interpolated by function values $u_{i, j+1}, u_{i+1, j+1}$ and $u_{i+2, j+1}$. Similarly, $u_{o, j}, u_{o, j-1}$ will be extrapolated by three function values from their right side in the $\Omega^{+}$ domain.

With the function values in $U_{o}$ appropriately approximated, (18) and (19) produce fictitious values $\hat{u}_{i, j}$ and $\hat{u}_{i+1, j}$ at points $(x_i, y_j)$ and $(x_{i+1}, y_j)$ represented as a linear combination of surrounding points near the interface and three jump conditions at interface point $(x^*, y^*)$. For instance, the representation of $\hat{u}_{i, j}$ can take a general form as

where $\mathbb{S}_{i, j}$ represents a set of surrounding nodes involved in (18) and (19).

In this manner, the needed two layers of fictitious values can be generated along $x-$direction by enforcing the jump conditions (16) or along $y-$direction by enforcing (17). Interested readers are referred to [42] for more details.

Remark 1. The above MIB formulation utilizes extrapolations in $\Omega^+$ only, i.e., $u_y^{+}$ and $u_x^{+}$ in (16) and (17) are perpendicularly approximated by using auxiliary values, and each auxiliary value will be extrapolated by using function values in $\Omega^+$. Recall that Eqs.(16) and (17) are derived after $u_y^-$ or $u_x^-$ being eliminated from Eqs.(13) - (15). On the other hand, if we remove $u_y^+$ or $u_x^+$ from Eqns (13) - (15), we can have another set of formulations involving $u_y^-$ or $u_x^-$. Numerically, $u_y^-$ or $u_x^-$ needs to be extrapolated by only using function values in $\Omega^-$. We call the first approach as MIB-PLUS and second one as MIB-MINUS in terms of function values used from $\Omega^{+}$ or $\Omega^{-}$ in the extrapolation process.

Remark 2. It is found in our experiments that the accuracy of derivative jump reconstruction in Eq. (11) could be affected by the fictitious value generated by MIB-PLUS or MIB-MINUS depending on $\beta^{+}$ and $\beta^{-}$. When $\beta^+>\beta^-$, the accuracy could be better if MIB-PLUS is used. Similarly, when $\beta^+<\beta^-$, the accuracy could be better if MIB-MINUS is adopted. Nevertheless, for efficiency concern, we will use fictitious values obtained by MIB-PLUS in the case when $\beta^{-}>\beta^{+}$, while MIB-MINUS for the case with $\beta^{-}<\beta^{+}$ [16]. This strategy can reduce the iteration number in the overall algorithm.

2.5.   Formulation of augmented system

Great details have been given to demonstrate how corrected differences are used to approximate the Laplacian operator. The right side of the Crank-Nicolson scheme (6) can be easily determined once the information of solution $u^n$, the interface conditions $\phi, \rho, \psi$ at current step $t^n$, and the source term $f$ at $t^n$ and $t^{n+1}$ are gathered. Especially, corrected differences with the jump approximation method (11) are utilized to deal with approximation at irregular points on both sides of the time-stepping scheme (6). However, it will be inefficient to apply iterative algorithms to solve the linear system underlying (6) at $t^{n+1}$. It is noted that the standard second order central difference matrix structure with some additional diagonal entries is obtained if the corrected difference is neglected. Due to such special matrix structure, a multigrid-based fast solver is established to fulfill the fast solution of the parabolic interface problem. Prior to the introduction of fast multigrid solver, we need to introduce the auxiliary variables in order to take advantage of the special matrix structure.

2.5.1.   Auxiliary variables

The right hand side of the time-stepping scheme (6) can be explicitly calculated with the known information on boundary or interface conditions, while the left hand side of (6) needs efficient algorithm design. For this purpose, we introduce auxiliary variables to form an augmented system.

With the general formulation of fictitious values in (20), the approximation (11) can be well defined. At each intersection point between the interface $\Gamma$ and $x-$or $y-$grid line, auxiliary variables $[u_{x}], [u_{xx}]$ or $[u_{y}], [u_{yy}]$ are introduced. Then approximation (11) can be reformulated as

where $C_k$ is the corresponding weight coefficient, $U$ signifies the needed function values in (11), and $\Phi_k$ is from the known quantities of the fictitious values in (11). The subscript $k$ indicates the jump of $k$th derivative at time $t^{n+1}$. Then the derivative jumps at all intersecting points can contribute a general matrix-vector form

where $C$ stands for the weight coefficients as in (21), $Q$ denotes all the auxiliary variables.

2.5.2.   Augmented system

Let $U_{i, j}$ indicate the discrete solution while $u(x_i, y_j)$ is continuous solution at $(x_i, y_j)$. Based on the corrected difference analysis for spatial discretization, the scheme (6) is expanded at all interior nodes as

where $C_{i, j}$ is the correction term, and $L_hU_{i, j}$ is the standard five points central difference scheme

Note that the subscript $(i, j)$ indicates the gird node $(x_i, y_j)$, and the superscript $n$ and $n+1$ imply the time $t^{n}$ and $t^{n+1}$, respectively. The degree of freedom of this system is $N = (n_x+1)(n_y+1)$. The correction term $C_{i, j}$ only exists at irregular points but vanishes at regular points. The discretization (23) gives rise to the matrix-vector form

where $D$ is a matrix with non-zero entries only along the diagonal, $\hat{A}$ is a $N$-by-$N$ non-symmetric and diagonally dominant matrix, the unknown $U$ is for solution at all grid nodes and $Q$ is for all the needed auxiliary variables. The diagonal entries of matrix $D$ are non-constant due to the different values of $\beta$ in different subdomain. As the right hand side of (26) can be determined in each time step, we can denote the right hand side of the above system as $F$. The above system can be further simplified as

after we combine the coefficient matrix $D$ and $\hat{A}$ into $A$, i.e., $A = D+\hat{A}$, and ignore the superscript $n+1$. Both $U$ and $Q$ are to be solved at time $t^{n+1}.$ Supposed $Q$ is already determined, then the time stepping process continues by updating right side $U^{n}$ with obtained solution $U^{n+1}$ in (26). The Dirichlet boundary condition is imposed in the vector $F$ with some rows of matrix $A$ composed by vector $(0, 0, \cdots, 0, 1, 0, \cdots, 0)$, i.e., only the diagonal element is one while other entries are zero. The total number of auxiliary variables $Q$, denoted as $M$ here, is double of that of all intersection points, while the latter is usually proportional to $n_x$ or $n_y$. In other words, $M$ is usually one-dimensional smaller than $N = (n_x+1)(n_y+1)$. The matrix $B$ with dimension $N$-by-$M$ contains coefficients from correction terms. In particular, matrix $B$ is a sparse matrix, since the correction terms only have impact on the irregular points which account for a small portion in the whole computation grid points.

We need to determine $Q$ first to solve for $U$ in (27). It is necessary to combine (27) and (22) to form augmented equation system

where

Note that the dimensions of matrix $C$ and $I$ match those of $A$ and $B$. Here $C$ is a $M$-by-$N$ sparse matrix, and $I$ is a $M$-by-$M$ identity matrix.

The matrix structure of $A$ allows for fast inversion with aid of multigrid method, which will be discussed later. Bearing this in mind, we will formulate two efficient solvers for the parabolic interface problem by using the following Schur complement strategy.

2.5.3.   Schur complement

In each time step, we aim at solving $U$ from

after we determine $Q$, and move all the known quantities to right hand side. Then the solution $U$ could be quickly obtained by applying the multigrid solver on the above equation thanks to the matrix structure of $A$. Such step is preceded by the determination of $Q$ from the Schur complement

which is obtained by eliminating $U$ from the augmented system (28).

Note that the linear system (30) for $Q$ has a much smaller degree of freedom than $U$, i.e. $M$. We propose two approaches to solve $Q$ from Eq. (30) in terms of iterative solver or the LU decomposition on the coefficient matrix of the left hand side in Eq.(30). The following steps are used to illustrate the procedure. One is by iterative algorithm of GMRES method(GMRES approach), while the other is by LU decomposition(LU approach).

● GMRES approach.

1. As a linear system, the right hand side RHS = $\varPhi-CA^{-1}F$ is determined by applying multigrid solver on $A^{-1}F$ and some arithmetic operations.

2. For the left hand side(LHS), we may consider an iteration approach by the GMRES method. In this approach, the matrix vector product of $(I-CA^{-1}B)Q$ is fulfilled in a couple of steps. Due to its equivalence to $I Q-CA^{-1}B Q$, a multigird solver is carried out on the product of $BQ$, i.e., $A^{-1}(B Q).$ Then it is followed by some more arithmetic operations on $I Q-CA^{-1}(B Q)$.

3. For the first time step $t = 0$, an initial guess $Q = (0, 0, \cdots, 0)^T$ is used to start the GMRES iteration. For the rest time steps, the initial guess of $Q$ is taken as the computed $Q$ from the previous time step. Either the maximal iteration number $5000$ or error tolerance equal to $10^{-10}$ terminates the GMRES iteration. The termination criteria are flexible according to actual computational needs.

● LU approach.

1. The first step is the same as that in the above GMRES approach.

2. For the left hand side(LHS), we will explicitly construct the matrix $I-CA^{-1}B$, where the inverse of A, i.e. $A^{-1}$, should not be formed, but is applied on the columns of matrix $B$. Then the explicit matrix is decomposed into the matrix product of $X$ and $Y$ from a LU decomposition solver. Such matrix construction and decomposition need to conducted only once before the time stepping.

3. The solution $Q$ in each time step can be simply obtained via the LU forward and backward substitutions.

Remark 3. The efficiency of GMRES approach depends on the iteration number for the GMRES iteration process. Note that each iteration involves one multigrid inversion with a complexity $O(N) = C_2N$ for some constant $C_2$. Fortunately, the iteration is for solving (30), whose dimension $M \ll N$. Generally, for any iteration solver, the better conditioned the iteration matrix is, the smaller the iteration number is. The selection of MIB-PLUS or MIB-MINUS mentioned in Remark 2 provides a good conditioning strategy. This makes the iteration number increase moderately as the mesh is refined.

Remark 4. The strategy of LU approach is different from GMRES approach in the sense that the Schur complement matrix is formed with multigrid algorithm applied for each column of matrix $B$. The more grid lines intersect the interface, the more columns matrix $B$ has. The computing cost is mainly contributed by the multigrid inversion on all the columns, before the time stepping. In the time-stepping process, only the forward and backward substitution of $LU$ are applied to solve for auxiliary variable, and one multigrid inversion is deployed on (29) for solution of $u$. Moreover, the LU approach is not impacted by the conditioner number of the linear system as it does for GMRES approach, since the LU decomposition is a direct rather than an iterative solver.

Remark 5. The complexities of the two approaches are in the same magnitude. As it is second order accurate for both temporal and spatial discretization, it is reasonable to set time increment $\Delta t$ proportional to spatial partition size $h$ to ensure the overall second order convergence. Suppose the total degree of freedom $N = n^2$, where $n$ is the grid number in each direction. The intersection number $M$ is hence equal to $C_0n$ for some constant $C_0$, and the temporal evolution number $N_t = C_3n$ for another constant $C_3$. We may roughly have an estimate of the complexities of the two approaches. The complexity of the LU approach is

where the first and second term in the first equality indicates the total sum of matrix generation cost arising from multigrid solver and the cost of multigrid for inversion deployed on (29) and right side of (30), and the LU forward and backward substitution cost in all time evolution steps. Besides, the last term in the first equality denotes the LU decomposition for the generated matrix in the first time stepping.

The computational complexity of the GMRES approach is

where $C_4$ indicates the sum of total times multigrid solver is call in the GMRES subroutine, and multigrid for inversion deployed on (29) and right side of (30).

The computational complexities of the two approaches are all of $O(n^3)$. Note that the number $C_4$ is a constant under the assumption that the iteration number of GMRES solver does not depend on mesh size. But in the practical implementation of our algorithm, the iteration number depends on the mesh size weakly. This can be seen in the numerical examples of the latter section. This makes AMIB-LU slightly more efficient than AMIB-GMRES when achieving second order accurate solution.

Remark 6. When an explicit finite difference method is employed to solve parabolic PDEs, a stability condition with $\Delta t < C h^2 / \beta_{max}$ has to be satisfied, where $\beta_{max} = \text{max}_{x\in \Omega}\beta(x)$. This means that one has to take $N_t = O(n^2)$ to ensure stability. It is reasonable to assume the complexity in each time step of this explicit scheme to be on the order of $O(n^2)$, since there are $n^2$ degree of freedom. Thus, the overall complexity of explicit finite difference solution of the 2D parabolic interface problems is about $O(n^4)$. This means that the proposed AMIB methods could be asymptotically faster than the explicit schemes.

2.5.4.   Multigrid algorithm

In the preceding section, a multigrid approach has been assumed in the discussion of building efficient algorithm. In this section, details on this solver are provided.

The multigrid approach is an efficient algebraic algorithm for solving a linear system

with known matrix $A$ and vector $B$. Recall in our study, $A$ has a special matrix structure: $A = D+\hat{A}$. Thus, Eq. (31) can be regarded as obtaining from a discretization to negative of Laplacian operator plus solution with variable coefficient as below:

where the domain and its partition are the same as that for the problem in this work, the discrete values of $a(x, y)$ on the given partition are the same as diagonal values of $D$, $U$ is the discrete solution to $u$ in (32), and function $b$ has its discrete data equal to $B$ at the interior grid nodes and for boundary conditions.

To be more clear, at the interior node,

for $1\leq i\leq n_x-1,$ and $1\leq j\leq n_y-1,$ while at the grid node on boundary of $\Omega$,

for $(i, j)$ such that $(x_i, y_j)\cap \partial \Omega \neq \emptyset.$

The discrete problem can then be reformulated as below according to the discrete data,

where $L_h$ and $D_h$ are the discrete versions of Laplacian operator, and the coefficient function $a(x, y)$ as defined in last section, and $B_h$ is the boundary operator. The $L_h$ and $D_h$ in (34) share the same stencil structure as in (24) and (25). When the mesh is refined by doubling $h$ to $H$, the value of $b_H$ will be sampled directly on the refined interior node, but the boundary value of $b_H$ is collected in the refined grid nodes on the boundary in the same way. The corresponding coefficient matrix for the discrete governing equation (34) on coarse mesh $H$ is a submatrix of that for mesh $h$.

This makes it possible for applying multigrid method to solve the linear system (31) efficiently. Once $U_h$ is determined, the solution of (31) is set as $U = U_h.$

The main idea of multigrid method lies in reducing the high frequency components of solution error on fine mesh, and that the low frequency errors are solved on coarse grids. Due to the fact that low frequency components of the errors can not be effectively reduced, the residual is restricted to coarser grids such that the low frequency errors are solved therein. Such determined errors are then interpolated onto fine meshes to provide some correction to the already obtained solution on fine meshes. Over the process, relaxation is needed on fine mesh to improve the solution accuracy, while restriction operator is imposed from fine mesh to coarse mesh and interpolation is needed from coarse mesh to fine mesh. Such process is iterated until the solution meets certain accuracy.

The relaxation is used to reduce errors in two aspects. One is for smoothing the high frequency component on the fine mesh. The other is used after the obtained errors are interpolated to add correction on solution on fine mesh. In this case, the relaxation is to reduce the error introduced from interpolation. In this work, we adopt the red-black Gauss-Seidal smoother as the relaxation scheme. The two-grid correction scheme as the simplest multigrid adopts only two level of meshes to achieve the solution correction via iteration between restriction and interpolation. We use such two-grid correction scheme to illustrate the way we apply multigrid scheme to obtain the solution of the linear system (31).

The above process continues until convergence. In our multigrid approach, we generalize the above process to more grid levels to avoid the exact correction solution for Eq.(36) in case that the cost is still very high on that coarse grid. This promotes the necessity of even coarser grids for correction solutions. The above procedure demonstrates the two-grid correction scheme using $V$-cycle. In our algorithm, we deploy $V$-cycle for the multigrid solver rather than the full multigrid approach since $V$-cycle tends to be more efficient in our practical implementation.

Restriction and interpolation are two needed transfer operators in the correction procedure. In our implementation, it is sufficient to adopt the following standard full-weighting stencil for restriction,

and the following standard interpolation scheme

The subscript $h$ and superscript $H$ in $R_h^H$ indicate restriction from fine mesh of mesh size $h$ to coarse mesh $H$, where $H = 2h.$ As shown in Fig. 4, the restricted function value at $(i, j)$ on coarse mesh comes from the nine function values at the highlighted nine grid points on fine mesh with the corresponding weights in matrix $R_h^H$. On the other hand, the subscript $H$ and superscript $h$ in $I_H^{h}$ indicate interpolation from coarse mesh $H$ to fine mesh of mesh size $h$, where $H = 2h.$ Through weight matrix $I_H^{h}$, function value at $(i, j)$ contributes to nine interpolated function values at the nine grid points shown in Fig. 4. To be clear, the value at $(i, j)$ on coarse mesh is multiplied respectively by the nine weights in matrix $I_H^{h}$ to give a contribution to the to-be-determined values at the nine grid points on fine mesh.

For the practical implementation of our multigrid algorithm for problem (34)(35), we set both $v_1$ and $v_2$ equal to 2. For the exact solution on the coarsest grid, we use biconjugate gradient iterative solver [31]. The iteration stopping criteria is set as error tolerance equal to $10^{-10}$, or maximal iteration number equal to $5000$. Besides, the multigrid stopping criteria is the residual $L_2$ error less than $10^{-10}$. We may use number of grid levels flexibly for the multigrid algorithm based on the mesh refinement.

3.   Numerical experiments

In this section, numerical experiments are carried out to demonstrate the accuracy, efficiency and stability of the proposed two AMIB algorithms for solving parabolic interface problems. For sake of notation, we name the AMIB method coupled with the GMRES approach on Schur complement as AMIB-GMRES, while the AMIB method associated with LU decomposition is called as AMIB-LU. For comparison purpose, the standard MIB method [42] with Crank-Nicolson scheme for time-stepping is also implemented.

For simplicity, a square domain with uniform number of partitions in each direction is assumed. We consider the number as $n_x = x_y$ in each direction. Time integration will be carried out from $t = 0$ until a stopping time $t = T$ with $n_t$ steps and $\Delta t = T/n_t$. The numerical accuracy and convergence will be tested by comparing the numerical solutions with exact solutions under the $L_{\infty}$ and $L_2$ norm.

The iteration solver used in MIB method is biconjugate gradient iterative solver [31]. The maximal iteration number is 5000, and the tolerance error is set as $10^{-10}.$ The number of mesh levels in the multigrid algorithm for AMIB is chosen as $4$ for $n_x = 16$, as $5$ for $n_x = 32$, as 6 for $n_x = 64$, and as 7 for all $n_x>64.$ For multigrid algorithm application, we have to set $n_x$ to be $2^k$ for a positive integer $k.$

All the experiments were carried out on MacBook Pro with 8.00 RAM and Intel 2.3 GHz Intel Core i5.

3.1.   Numerical accuracy

Example 1. In this example, we consider a parabolic problem on a square domain $[-\frac{\pi}{3}, \frac{\pi}{3}]\times[-\frac{\pi}{3}, \frac{\pi}{3}]$ with a circular interface defined by $r^2: = x^2+y^2 = \frac{1}{4}$. The exact solution to this problem is

with the diffusion coefficients

Here we call $\beta$ as $\beta^{+}$ when it is outside of $\Gamma$, and $\beta^{-}$ when it is inside of the interface $\Gamma.$ The interface jump conditions on an interface point $(x, y) = (\frac{1}{2}\cos(\theta), \frac{1}{2}\sin(\theta))$ can be derived as below:

It follows the same principle in below notations. The source term $f(x, y)$ is related to the above designated solution,

The initial condition and boundary conditions can be easily determined by the analytical solutions. Beside, we may observe that both the zeroth and first order jump conditions on the interface are time-independent. To test spatial convergence rate, we consider the fixed stopping time at $T = 1$ with the time increment $\Delta t = 10^{-3}$ being small enough such that temporal error will not interfere with the examination of the spatial convergence. Table 1 clearly reveals that the spatial discretization reaches the second order convergence in space.

On the other hand, the temporal convergence needs to be tested. For such concern, the spatial partition in each direction is set to be $n_x = 512$, and various time steps are deployed. The results in Table 2 indicates that the second order temporal convergence can be monitored from the second row to the last second row. Such second order convergence starts to occur from $\Delta t = 0.25$. The total CPU time of each time evolution methods is listed in the last column in Table 2, where we can compare the efficiency between these method. It can be observed that if the number of time steps is small, AMIB-GMRES is superior than AMIB-LU. However, in the long run, AMIB-LU outperforms both AMIB-GMRES and MIB methods. Take the case of $n_t = 128$ for instance, CPU time of MIB is about 12 times larger than that of AMIB-LU. Similar results can also be found in the last row of Table 1 of spatial convergence test when partition number is $512$ in each direction. The number of time evolution steps is 1000, leads to CPU cost of MIB being 14 times larger than that of AMIB-LU. The inefficiency of MIB or other interface algorithms implemented with iterative solvers for parabolic interface problems leads to the demand for the development of efficient solvers.

Example 2. In this example, we consider a parabolic problem on a square domain $[-\frac{\pi}{3}, \frac{\pi}{3}]\times[-\frac{\pi}{3}, \frac{\pi}{3}]$ with a circular interface defined by $r^2: = x^2+y^2 = \frac{1}{4}$. The exact solution to this problem is

The source terms can be correspondingly determined,

The zeroth and first order interface can be found to be time and space dependent, representing the most general jump conditions. In this example, we consider $k = 2$.

It is also reasonable to approximate the equation (4) with

with the difference to Eq.(5) on the approximation to $f$ at time $t = t_n+0.5\Delta t$ by $f_{i, j}^{n+\frac{1}{2}}$ in Eq.(41), while $\frac{1}{2}(f_{i, j}^{n}+f_{i, j}^{n+1})$ is used in Eq.(6). We are interested in the numerical impacts with such different treatment for the source term $f.$ Table 3 shows their numerical comparisons between scheme (41) and (5) with $f^{n+\frac{1}{2}}$ indicating scheme (41) and the other part for scheme (5). In this experiment, the partition number in each direction is $n_x = 512.$ The temporal convergence can all be determined according to the results. It can be found that the solution accuracy from (5) with $n_t = 128$ is comparable to that of (41) with $n_t = 256$. Such accuracy improvement justifies that scheme (5) is superior to (41) in terms of accuracy when dealing with the interface jump conditions that are time and space independent. On the other hand, it can be found in Tables 4 and 5 that the time costs from scheme (5) and (41) are basically the same for each $n_t$. From this point of view, scheme (5) is a better choice for time evolution than (41). Thus we will focus on scheme (5) in the following discussion to deal with the parabolic interface problem under all kinds of circumstances.

The error versus time partition increment $\Delta t$ on log-log scale from scheme (5) is plotted on Fig. 5, from which we can observe that the temporal second order accuracy is obtained for both AMIB-GMRES and AMIB-LU methods.

On the other hand, the spatial convergence are investigated on the two methods. The stopping time is fixed as $T = 2$, and the time increment is set as $0.002.$ The second order spatial convergence of the two approaches concerned with source term $f$ can be observed from the results in Table 4. The comparison shows no much difference on spatial accuracy by scheme (5) and (41).

We then want to study the impact of ratio of $\beta^+:\beta^-$ on the stability for scheme (5). We will fix one of two $\beta$ values equal to 1, then vary the value of the other one to make the ratio large or small enough. We fix the partition number in each direction as $128$, and the stopping time equal to $2$ with time increment $0.002$. It can be observed from the Table 5 that the two AMIB methods are all stable enough to deal with large contrast of the coefficients. On the other hand, the varying coefficient contrasts seem to have larger impact on the CPU time of AMIB-GMRES than AMIB-LU. The CPU time of AMIB-GMRES become twice expansive when the ratio of $\beta^+:\beta^-$ is raised from $10:1$ to $10000:1$. This is also the same case for AMIB-GMRES when ratio of $\beta^+:\beta^-$ is raised from $1:10$ to $1:10000.$ The milder impact on AMIB-LU is due to the direct solver by LU decomposition rather than iterative solvers whose conditional number can be affected by the ratio of $\beta^+:\beta^-$. Such conclusion is based on the comparison between the two AMIB methods vertically in Table 5.

Note that the "Small $\beta$ side" indicates the the MIB extrapolation part is fulfilled by the function values from the domain related to the small $\beta.$ The "Large $\beta$" side indicates the domain related to the large $\beta$ in a similar fashion. From the horizontal comparison for each AMIB method in Table 5, MIB fictitious values with the "Small $\beta$ side" extrapolation is superior to the one with the "Large $\beta$ side", because the CPU time is less than that in the other way around for AMIB-GMRES. This justifies the selection of $\beta$ in remark 2 for concern of being well-conditioned. Hence, it is preferable to adopt the fictitious values with the extrapolation part fulfilled using the function values from domain of "small $\beta$." The accuracy of the approximate solutions concerned with the "Small" or "Large" $\beta$ has no much big difference using AMIB-GMRES or AMIB-LU in this example. In the following examples, we will adopt the "Small $\beta$ side" strategy for both AMIB schemes.

Example 3. We are interested in the performance of the two AMIB methods on complicated interface shape. In this example, we focus on the case of interface shape governed by parametric function

embedded in the domain $[-\frac{\pi}{3}, \frac{\pi}{3}]\times[-\frac{\pi}{3}, \frac{\pi}{3}]$. The solution and all the related boundary and interface conditions are set as in the last example. The stopping time is given as $T = 1$, and the coefficients are $(\beta^+, \beta^-) = (100, 1).$ The spatial accuracy and convergence are examined with the small time increment $\Delta t = 0.001.$ The numerical results in Table 6 on the successively refined mesh demonstrate the second order spatial accuracy.

We next examine the temporal order of the two AMIB methods by setting $n_x = 512$, and the stopping time as $T = 1.$ As the time evolution scheme (5) is efficient enough for the temporal convergence, the numerical results in Table 7 show that both AMIB methods obtain their limiting convergence very quickly, leading to the temporal second order observed from $n_t = 2$ to $n_t = 4$, and the reduced order in further steps.

Since the AMIB algorithm achieves second order accuracy for both spatial and temporal discretization, we test the computing efficiency by refining the spatial and temporal steps simultaneously. See the Table 8 with stopping time $T = 1.$ We can observe that the second overall accuracy of the three methods are all achieved after the spatial and temporal partition are refined at the same time. Such detected second order accuracy can validate the temporal second order despite the reduced order observation in the Table 7. The CPU time comparison of the three methods further demonstrate the better efficiency of the two AMIB method over MIB method. Besides, the time ratios of MIB method with mesh refinement are higher than the two AMIB methods, meaning the MIB is the most expansive algorithm among the three. AMIB-LU seems a little more efficient than AMIB-GMRES thanks to its relatively smaller time ratio when the mesh is refined. The reason for its better efficiency is mainly due to a smaller average number of multigrid solver (ANMG) being used in each time step as recorded in the last column in Table 8. It can be seen that the average number for AMIB-GMRES mildly increases as the mesh is refined. This is due to the fact that the iteration number of GMRES solver in solving (30) weakly depends on the grid numbers. However, the average number remains constant for AMIB-LU, since the majority of the number of multigrid solver used is on forming the explicit matrix of left side of (30).

It can be noticed that the numerical results for MIB in the last row of Table 8 are not available due to failure of convergence from the biconjugate gradient solver [31]. The non-convergence is not caused by stability issue on the algorithm, as a different iterative solver, such as GMRES, can produce convergent iterative solution in our numerical test. This demonstrated that the performance of the AMIB methods does not heavily depend on the iterative solvers, while the MIB method does.

Example 4. We are interested in the stability performance of the two methods on complicated interface shape. In this example, we focus on the case of interface shape governed by parametric function

embedded in the domain $[-\frac{\pi}{3}, \frac{\pi}{3}]\times[-\frac{\pi}{3}, \frac{\pi}{3}]$.

The exact solution to this problem is

with the diffusion coefficients

The spatial convergence is validated by refining the mesh with the stopping time $T = 2$, and time increment $\Delta t = 0.002.$ In Table 9, the second order accuracy can be observed for both AMIB methods.

We next examine the temporal order of the two AMIB methods by setting $n_x = 512$, and the stopping time as $T = 0.001.$ The second order temporal convergence can be observed in Table 10 from $n_t$ equal to 2 up to 16.

To numerically analyze the stability of the two AMIB schemes, we examine different time stepping. We choose $n_x = 128$ and stopping time $T$ equal to $10^4\Delta t$ as shown in Fig. 6. All numerical solutions are found to be stable with large temporal increments. And the corresponding numerical errors are bounded as shown in Fig. 6. This justifies the reliable stability of our methods for solving parabolic interface problem.

Example 5. We next study a problem with multiple subdomains. Consider a 2D parabolic equation

defined in a square domain $\Omega = [-\frac{\pi}{3}, \frac{\pi}{3}]\times[-\frac{\pi}{3}, \frac{\pi}{3}]$, which is divided into three regions by two ellipse shape interface:

Three pieces of solutions are defined in each subdomain

where parameter $k$ is set to be 5. The diffusion coefficients are also defined respectively as

and the parameter $k$ is set as 5. The source term $q(x, y)$ can be determined by the analytical solutions, and Dirichlet boundary conditions are assumed.

The second order spatial and temporal accuracy of solutions have been examined through above examples. Besides, we want to further investigate the convergence rate of the gradient recovery. Second order central difference is adopted to approximate the gradients with fictitious values replacing obtained numerical real function values when the stencil cuts through the interface. Table 11 shows that second order convergence is achieved for both solutions and gradients from the three methods in solving the multi-domain problem. With the confirmed second orders of the three methods, we pay attention to the comparison of the methods on overall efficiency for solving multi-domain problem. Table 11 demonstrates the efficiency test, where it can be observed that AMIB-LU outperform the other methods for solving such problems. The last column of Table 11 shows that AMIB-LU has smaller time ratio than AMIB-GMRES or AMIB with mesh refinement. In particular, the computing time by MIB is almost 19 times larger than that of AMIB-LU in dealing the multi-domain problem on the $(n_x, n_t) = (513, 160)$ mesh, further demonstrating the efficiency improvement from AMIB-LU.

We note that there exist more challenging multi-material interface problems, in which two or more material interfaces join together or cross each other. A second order accurate MIB method has been developed in [38] for elliptic equations with intersecting interfaces. Based on such a MIB method, an augmented MIB method could be developed for intersecting interface problems. This is an interesting direction to explore in the future.

Example 6. We investigate a physical problem modeled by the heat equation without analytical solutions

subject to the homogeneous interface condition

on the interface $r^2: = x^2+y^2 = \frac{1}{4}$ embedded in the domain $[-\frac{\pi}{3}, \frac{\pi}{3}]\times[-\frac{\pi}{3}, \frac{\pi}{3}]$. The diffusion coefficients

The Dirichlet boundary condition equal to zero is imposed at the boundary of given domain. Besides, the initial condition is defined as

with parameter $\sigma = 0.1.$ As the analytical solution is not readily available, we generate the reference solution on certain mesh to examine the numerical performance of the two AMIB methods.

The spatial convergence is validated by refining the mesh under time increment $\Delta t = 0.001$ for stopping time $T = 1.$ The reference solution is created on the mesh $(n_x, \Delta t) = (512, 0.0005)$. The second order accuracy can be observed for both AMIB methods in Table 12.

The temporal convergence is validated in a similar fashion by refining time increment for a fixed stopping time $T = 1$, and the spatial discretization is based on mesh $(n_x, n_y) = (256, 256).$ The reference solution is obtained on the mesh $(n_x, \Delta t) = (512, 0.0005)$. Due to the practical physical phenomenon modeled by heat equation, the temperature quickly drops in a very short time, leading to large approximation errors when a large time step size is adopted as shown in Fig. 7. Nevertheless, small time steps can produce well approximated solutions. This is true for both AMIB-GMRES and AMIB-LU methods, especially when $n_t>160$ is employed. For such reason, it is difficult to observe the expected second order temporal convergence.

The figure 8 shows solution evolution over time for this physical model. As discussed before, the temperature drops quickly in a short time. This can be observed in Fig. 7 when time is changing from 0 to 0.1. When time equals to 1, the temperature has dropped to a magnitude of $10^{-4}.$

4.   Summary

In this work, a new finite difference method has been proposed to solve parabolic interface problem with second order accuracy achieved in both time and space. The efficiency is fulfilled with multigrid method built in the Schur complement process in the framework of an augmented system. A simple multigrid solver is designed to solve the linear system due to the special structure of coefficient matrix as part of the augmented system. Two approaches based on the multigrid solver are proposed in the time-step process, with one incorporated in GMRES iteration for each time evolution, and the other one used in the coefficient matrix formation of LU decomposition solver in the first time step. The latter approach with LU decomposition can be more efficient especially when a lot of time-stepping evolutions are needed for solving the given parabolic interface problem. Both AMIB methods are significantly faster than the standard MIB method. Moreover, being free of stability constraint, the present implicit method could be asymptotically faster than explicit schemes.

The immersed interface method (IIM) has been applied to variable coefficient elliptic equations [4,21], and these studies demonstrated how to adopt numerical components, such as jump-corrected differences, augmented system, and multigrid solver, for variable coefficient problems. Thus, based on these studies, it is possible to generalize the present mulitgrid AMIB method for variable coefficient interface problems. This will be explored in the future.

Acknowledgments

This research is supported in part by the National Science Foundation (NSF) of USA under grant DMS-1812930.