A high order immersed finite element method for parabolic interface problems

We present a high order immersed finite element (IFE) method for solving 1D parabolic interface problems. These methods allow the solution mesh to be independent of the interface. Time marching schemes including Backward-Euler and Crank-Nicolson methods are implemented to fully discretize the system. Numerical examples are provided to test the performance of our numerical schemes.


Introduction
We consider the following one dimensional parabolic interface problem: u(t, x) = g(t, x), x ∈ ∂Ω, t ∈ (0, T ], u(0, x) = u 0 (x), x ∈ Ω. (1) Assume that Ω is an open interval separated by an interface point α. The interface point divides Ω into two sub-domains Ω + and Ω − such that Ω = Ω + ∪ Ω − ∪ {α}. We assume that there is only one type of material in each sub-domain. This means that the coefficient function β is continuous within each sub-domain but may be discontinuous across the interface: The solution u is assumed to satisfy the following interface jump conditions: For the sake of simplicity, we consider only one interface point. The extension is straightforward for the case of multiple interface points.
The interface problems (1)-(4) arise in many applications in science and engineering. A wide variety of numerical methods have been developed to solve PDE interface problems, such as the finite element method [1]. However, the convergence rate may not be optimal unless the interface is resolved by the solution mesh. In other words, the mesh needs to fit the interface. For problems that involve a moving interface, the body fitting restriction requires a new mesh generation at each time level. This will be time consuming, especially for the complicated interface shapes in the 2D or the 3D domain. Consequently, the overall computational cost will increase.
Research has been done to develop numerical methods based on interface independent mesh, such as immersed finite element methods (IFEM) [2,3]. The basic idea is to modify the finite element basis functions locally to satisfy the jump conditions associated with the interface. This method has been extended to solve moving interface problems [4][5][6]. Most of the existing literature for time-dependent problems are for lowest order space approximations. In this paper, we consider some higher-order IFE approximations for one-dimensional parabolic interface problems. For temporal discretization, we utilize the standard Backward-Euler and Crank-Nicolson methods. Our numerical schemes converge in optimal orders in both spatial and temporal discretization even if our solution meshes do not fit the interface.
The rest of the article will be organized as follows. In section 2, we recall some preliminary results of the immersed finite element method. In section 3, we derive the semi-discrete and the fully-discrete numerical schemes for solving the parabolic interface problem. Lastly, in section 4, numerical examples are provided to demonstrate features of our schemes.

Preliminaries of the Immersed Finite Element Method
In this section, we briefly recall the one-dimensional IFEM spaces. The main idea of IFEM is to locally modify the FEM basis function to accommodate the flux jump condition (4). To be more specific, let us consider an arbitrary partition a = x 0 < x 1 < x 2 < · · · < x L = b of the domain Ω. For an interface problem, there exists at least one element containing the interface point. Without loss of generality, let us assume that the interface point α lies in the element K = (x 1 , x 2 ). This element is called an interface element, and the rest of elements are called non-interface elements. The local linear IFE space consists of two piecewise linear functions φ 0,K and φ 1,K . They are formed to satisfy the nodal value conditions and the interface jump conditions as follows For higher order IFE basis functions, we refer to [7] for Lagrange type IFE basis functions and [8] for orthogonal IFE basis functions. Plots for some orthogonal high order IFE basis functions are shown in Figure 1.
Following standard procedure, we can define the global IFE basis functions by requiring the continuity across the element boundaries. Let φ j , j = 1, · · · , L, be the global IFE basis functions. Then the global IFE space can be written as

Numerical Schemes
In this section, we present the IFE method of lines (MOL) algorithm for solving parabolic interface problems. The method of lines technique is useful since the spatial and the temporal discretizations are done separately. This allows MOL to be used in a variety of problems to achieve accurate numerical approximations.

Semi-Discretization
We first consider the spatial discretization for model problem (1). To derive its weak formulation, we multiply (1) by any test function v(x) ∈ H 1 0 (Ω). Applying integration by parts and the flux jump conditions, we have Note that the interface jump condition (4) must hold to obtain the weak formulation of the interface problem. Define the IFE trial and test function spaces as The semi-discrete problem for each time t is: find u h ∈ U h such that where This can be equivalently written in the matrix-vector form The notations in (9) are specified as follows • M is the mass matrix and • u(t) = [u 1 (t), u 2 (t), · · · , u N (t)] T is the unknown vector.

Full-Discretization
We further discretize the ODE system (9) by some well-known difference schemes such as Forward/Backward-Euler methods and Crack-Nicolson methods. We derive a general framework to incorporate all these methods by introducing a parameter θ ∈ [0, 1].
The Forward-Euler method is notoriously unstable as it creates heavy oscillations in general. On the other hand, both the Backward-Euler (BE) method and the Crank-Nicolson (CN) method are very stable. From the derivation, we can expect the convergence rates in the temporal discretization to be O(∆t) and O(∆t 2 ) for BE and CN, respectively. Combining with the optimal convergence for spatial discretization from Section 3.1, we look to achieve the overall convergence rate listed in Table 1. In this table, k denotes the polynomial degree for IFEM spaces.

Numerical Results
In this section, we present a few numerical examples to test the performance of our schemes.
For each test, we use both linear and quadratic IFEMs for spatial discretization on a family of uniform meshes {T h }, where each mesh size is h = 1 L . The temporal discretization uses the Backward-Euler and Crank-Nicolson methods on a uniform time partition with step size ∆t equal to the spatial mesh size h.
One can easily verify that the solution (14) satisfies the jump conditions (3) and (4). We report the error of numerical solution at the final time level in L ∞ , L 2 , and H 1 norms. Numerical results for linear and quadratic IFEMs with the Backward-Euler scheme are reported in Table 2 and Table 3, respectively. Numerical results for linear and quadratic IFEM with the Crank-Nicolson scheme are reported in Table 4 and Table 5, respectively. We observe from these tables that the convergence rates in all three norms match the expected convergence rates in Table 1.

Example 2: Parabolic Interface Problem With Non-Constant Coefficient β
In this example, we consider the interface problem with a large contrast non-constant coefficient function as defined below:  where the interface α = 5/6 in this example. The exact solution is Tables 6 and 7 show the error and convergence rates for the Backward-Euler solution with linear and quadratic IFEM functions, respectively. Tables 8 and 9 report the errors for Crank-Nicolson method using linear and quadratic IFEM functions. Again, these results are consistent with the anticipated convergence order in Table 1. This example shows the robustness of our numerical schemes with respect to the non-constant coefficient functions and high-jump circumstances.