# Stability and convergence analysis of a domain decomposition FE/FD method for the Maxwell's equations in time domain

Stability and convergence analysis for the domain decomposition finite element/finite difference (FE/FD) method is presented. The analysis is designed for semi-discrete finite element scheme for the time-dependent Maxwell's equations. The explicit finite element schemes in different settings of the spatial domain are constructed and domain decomposition algorithm is formulated. Several numerical examples validate convergence rates obtained in the theoretical studies.

## Authors

• 1 publication
• 2 publications
09/23/2020

### Finite element analysis for a diffusion equation on a harmonically evolving domain

We study convergence of the evolving finite element semi-discretization ...
11/24/2021

### Sensitivity Analysis of Domain Decomposition method for 4D Variational Data Assimilation (DD-4DVAR DA)

We prove consistence, convergence and stability of the Domain Decomposit...
04/12/2021

### Analysis of algebraic flux correction for semi-discrete advection problems

We present stability and error analysis for algebraic flux correction sc...
08/08/2017

### Modelling of a Permanent Magnet Synchronous Machine Using Isogeometric Analysis

Isogeometric analysis (IGA) is used to simulate a permanent magnet synch...
05/25/2021

### Parallel domain decomposition solvers for the time harmonic Maxwell equations

The time harmonic Maxwell equations are of current interest in computati...
05/15/2021

### A FETI approach to domain decomposition for meshfree discretizations of nonlocal problems

We propose a domain decomposition method for the efficient simulation of...
05/06/2020

### Optimally Convergent Mixed Finite Element Methods for the Stochastic Stokes Equations

We propose some new mixed finite element methods for the time dependent ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

New computational techniques meet the needs of industry in developing efficient computational methods to simulate partial differential equations (PDEs). Especially, for simulations in higher dimensions and large computational domains. In this regard, certain

domain decomposition (DD) method leading to efficient schemes, in numerical investigations, has gained a lot of interest in numerical analysis community. This variant of the DD method is the subject of our current, and some related and ongoing, research. This type of schemes was previously studied, e.g. in [12, 30] for others than studied here problems.

The present work is a further development of the DD hybrid finite element/finite difference (FE/FD) method for time-dependent Maxwell’s equations for electric field in non-conductive media, studied in [3, 4]. A stable, time-domain (TD), DDM scheme for Maxwell’s equations was proposed in [26, 27], and further verified in [13]. This method uses FDTD scheme on the, structured, FD part of the mesh, and edge elements on unstructured part. In applications, because of edge elements implementation, this method remains computationally expensive.

The domain decomposition FE/FD method for time-dependent Maxwell’s equations for electric field, assuming constant dielectric permittivity function in a finite difference domain, was considered in [3]

. This assumption simplifies the numerical schemes in both FE and FD domains and significantly reduces computational efforts for implementation of the whole DD method. Modified numerical scheme, energy estimate and numerical verifications of this method was presented in

[4]. However, the fully stability and convergence analysis with numerical implementations, in - and -norms of the developed FE and FD schemes, are not presented in the above studies. We fill this gap in the present work.

More specificaly, we present stability analysis for explicit schemes for both FEM and FDM in the DD hybrid FE/FD method. The DDM is constructed such that FEM and FDM coincide on the common, structured, overlapping layer between the two subdomains. The resulting domain decomposition approach at the overlapping layers can be viewed as a FE scheme which avoids instabilities at the interfaces. Similar to the DD approach of [4, 5], we decompose the computational domain such that FEM and FDM are used in different subdomains: FDM in simple geometry and FE in the subdomain where more detailed information is needed about the structure of this subdomain. This also allows application of adaptive FEM in such subdomain, see, e.g. [2, 3, 8, 9, 10, 28, 29].

Reliability and convergence of the domain decomposition method, studied in this work, are evident for solution of coefficient inverse problems (CIPs) in , see, e.g. [8, 9, 10, 28, 29]. For the case of CIPs, the computational domain is splitted into subdomains such that a simple discretization scheme can be used in a large region and more refined discretization scheme is applied in smaller, however more critical, part of the domain. In most algorithms for solution of electromagnetic CIPs, to determine the dielectric permittivity function inside a computational domain, a qualitative collection of experimental measurements is necessary on it’s boundary or in it’s neighborhood. In such cases it is convenient to condsider the numerical solution of time-dependent Maxwell’s equations in different subdomains with constant dielectric permittivity function in some subdomain and non-constant in the other ones. For the time-dependent Maxwell’s equations, the DD scheme of [4], which is analyzed in the present work, is used for solution of different CIPs to determine the dielectric permittivity function in non-conductive media using simulated and experimentally generated data, see [2, 8, 9, 10, 28, 29].

An outline of this paper is as follows. In Section 2 we introduce the mathematical model. In Section 3 we briefly present the domain decomposition FE/FD method and communication scheme between two methods. In Section 4 we describe the domain decomposition FE/FD method for solution of Maxwell’s equations and set up the finite element and finite difference schemes. Section 5 is devoted to the stability analysis. In Section 6, we derive optimal a priori error estimates in finite element method for the semi-discrete (spatial discretizations) problems. Finally, in Section 7 we present numerical implementations that justify the theoretical investigations of the paper. In what follows, will be a generic constant independent of all parameters, unless otherwise specifically specified, and not necessarliy the same at each occurance.

## 2 The mathematical model

The Cauchy problem for the electric field , , , of the Maxwell’s equations, under the assumptions that the dimensionless relative magnetic permeability of the medium is and the electric volume charges are zero, is given by

 1c2εr(x)∂2E∂t2+∇×∇×E=−μ0σ(x)∂E∂t,∇⋅(εE)=0, (1)

where, and are the dimensionless relative dielectric permittivity and electric conductivity functions, respectively. , and are the permittivity and permeability of the free space, respectively, and is the speed of light in free space. In this paper we consider the problem (1) in non-conductive media, i.e. , and hence study the initial value problem

 ε(x)∂2E∂t2+∇×∇×E=0,x∈R3,t∈(0,T].∇⋅(εE)=0,E(x,0)=f0(x),∂E∂t(x,0)=f1(x),  x∈R3,t∈(0,T]. (2)

To solve the problem (2) numerically, we consider it in a bounded domain (instead of whole ), with boundary , and employ a split scheme on : a hybrid, finite element/finite difference scheme, kind of domain decomposition, developed in [3, 4] and summarized in Algorithm 1. More specifically, we divide the computational domain into two subregions, and such that and is a subset of the convex hul of . The function is assumed to be constant in , and bounded and smooth in . The communication ß between and is arranged using an overlapping mesh structure through a two-element thick layer around as shown by blue and green common boundaries in Figure 1. The blue boundary is outer boundary of and inner boundary of . Similarly, the green boundary is the inner boundary of from which the solution is copied to the green boundary of .

The key idea with such a decomposition is to be able apply different numerical methods in different computational domains. For the numerical solution of (2) in we use the finite difference method on a structured mesh. In , we use finite elements on a sequence of unstructured meshes , with elements consisting of triangles in and tetrahedra in , both satisfying minimal angle condition. This approach combines the flexibility of the finite elements and the efficiency of the finite differences in terms of speed and memory usage and fits well for reconstruction algorithms presented below.

We assume that for some known constant , the function satisfies

 ε(x)∈[1,d],  for x∈ΩFEM,  ε(x)=1,  for x∈ΩFDM. (3)

Conditions (3) on and the relation

 ∇×∇×E=∇(∇⋅E)−∇⋅(∇E), (4)

together with divergence free field , make equations in (2) independent of each others in and so that, in , we just need to solve the system of wave equations:

 ∂2E∂t2−ΔE=0,   (x,t)∈Ω2×(0,T]. (5)

Remark

It is well known that, for stable implementation of the finite element solution of Maxwell’s equations, divergence-free edge elements are the most satisfactory ones from a theoretical point of view [21, 24]. However, the edge elements are less attractive for solution of time-dependent problems, since a linear system of equations should be solved at each time iteration step. In contrary, P1 elements can be efficiently used in a fully explicit finite element scheme with lumped mass matrix [11, 14, 18]. It is also well known that numerical solution of Maxwell’s equations using nodal finite elements is often unstable and results spurious oscillatory solutions [22, 25]. There are a number of techniques to overcome such instabilities, see, e.g. [15, 16, 17, 23, 25].

In [6, 7], a finite element analysis shows stability and consistency of the stabilized finite element method for the solution of (1) with . In the current study we show stability and convergence for the combined FEM/FDM scheme, under the condition (3) on , where the stabilized FEM is used for the numerical solution of (2) in and usual FDM discretization of (5) is applied in .

Remark

Here, we consider the case when for . Further, we assume that , and . Recall that we assumed non-conductive media: . In the presence of electric conductivity, additional -terms appear in the equations. They lead to more involved estimates and heavier implementations which we plan to perform in a forthcoming study.

Hence, in this note we study the following initial boundary value problem:

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩ε∂ttE+∇×∇×E=0 in Ω×(0,T),E(⋅,0)=f0(⋅), and ∂tE(⋅,0)=f1(⋅)% in Ω,E=0 on ∂Ω×(0,T),∇⋅(εE)=0 in Ω. (6)

## 3 The structure of domain decomposition

We now describe the DD method between two domains, and , where the FEM is used for computation of the solution in , and FDM is used in . Communication between and is achieved letting overlapping of both meshes across a two-element thick layer around - see Figure 1. The common nodes of both and domains belong to either of the following boundaries (see Figure 1):

• Nodes on the blue boundary - lie on the boundary of and are interior to ,

• Nodes on the green boundary - lie on the inner boundary of and are interior to .

Then the main loop in time for the explicit hybrid FEM/FDM scheme, that solves (2) associates with appropriate boundary conditions, at each time step is described in Algorithm 1 below:

By (3), at the overlapping nodes between and . Thus, FEM and FDM schemes coincide on the common, structured, overlapping layer. Hence, we avoid instabilities at interfaces.

## 4 Derivation of computational schemes

In this section we construct, combined, finite element-finite difference schemes to solve the model problem (6). To do this, first we present the finite element scheme to solve (6) in entire . This induces the finite element scheme in . Then, we derive the finite difference scheme in , when the domain decomposition FE/FD structure is applied to solve (6) in .

Remark

The computational schemes derived in this section are explicit and therefore, for their converegence, the CFL condition below (see ,e.g. [6, 7]) should be satisfied

 τ≤hη,   η=C√1+3∥ε−1∥∞. (7)

Here is a mesh independent constant, is the time step, and is the mesh size.

In the sequel, we denote the inner product of by , and the corresponding norm by . The scalar inner product in we denote by , and the associated norm by . Further, we let to be the boundary of , the inner boundary of and the outer boundary of .

### 4.1 Finite element discretization in Ω

First we derive finite element scheme to solve the model problem (6) in whole . Next, we discretize in two steps: (i) the spatial discretization using a partition of into elements , where is a mesh function defined as , representing the local diameter of elements. We also denote by a partition of the boundary into boundaries of the elements such that, at least one of the vertices of these elements belong to . (ii) As for temporal discretization, we let be a uniform partition of the time interval into subintervals of length As usual, we also assume a minimal angle condition on elements in . To formulate the finite element method for (6) in , we introduce the finite element space for each component of the electric field defined by

 WEh(Ω):={w∈H1(Ω):w|K∈P1(K),∀K∈Kh},

where denote the set of piecewise-linear functions on . Setting we define (resp. ) to be the usual

-interpolant of

(resp. ) in (6). Also, because of the Dirichlet boundary data in (6) we need to choose the test function space as

 WEh,0:={v∈WEh|v=0,on ∂Ω}. (8)

Then, recalling (4), and the fact that the spatial semi-discrete problem in reads:

Find such that ,

 (ε∂ttEh,v)+(∇Eh,∇v)+(∇⋅(εEh),∇⋅v)−(∇⋅Eh,∇⋅v)−⟨∂nEh,v⟩∂Ω:=∑5j=1Tj=0,Eh(⋅,0)=f0h(⋅) and ∂tEh(⋅,0)=f1h(⋅) in Ω. (9)

We note that

 T3=(∇⋅(εEh),∇⋅v)=(∇ε⋅Eh,∇⋅v)+(ε∇⋅Eh,∇⋅v),

implies that

 T3+T4=(∇⋅(εEh),∇⋅v)−(∇⋅Eh,∇⋅v)=(∇ε⋅Eh,∇⋅v)+((ε−1)∇⋅Eh,∇⋅v). (10)

Recalling (8), vanishing test functions at the boundary yields , and the final weak formulation for the semi-discrete problem in is: Find such that ,

 BΩ(Eh,v):=(ε∂ttEh,v% )+(∇Eh,∇v)+(∇ε⋅Eh,∇⋅v)+((ε−1)∇⋅Eh,∇⋅%v)=0. (11)

For the, reflexive, inhomogeneous boundary condition, see the FE scheme for .

To get fully discrete scheme for (6) we apply time discretization to (9) approximating , denoted by , where we use the central difference scheme, for :

 (εEk+1h−2Ekh+Ek−1hτ2,v)+(∇Ekh,∇v)+(∇⋅(εEkh),∇⋅v)−(∇⋅Ekh,∇⋅v)=0∀v∈WEh0(Ω),E0h=f0h and E1h=E0h+τf1h in Ω. (12)

Multiplying both sides of (12) by we get, that

 (Ek+1h−2Ekh+Ek−1h,v)+τ2(1/ε∇Ekh,∇v)+τ2(1/ε∇⋅(εEkh),∇⋅v)−τ2(1/ε∇⋅Ekh,∇⋅v)=0,E0h=f0h and E1h=E0h+τf1h, in Ω. (13)

Rearranging terms in (13) we get the following scheme: Given the approximate initial data, and , find such that ,

 (Ek+1h,v)=(2Ekh,v)−(Ek−1h,v)−τ2(1/ε∇Ekh,∇v)−τ2(1/ε∇⋅(εEkh),∇⋅v)+τ2(1/ε∇⋅Ekh,∇⋅%v)E0h=f0h and E1h=E0h+τf1h in Ω. (14)

### 4.2 Finite element discretization in ΩFEM

To solve the model problem (6) via the domain decomposition FE/FD method, we use the split , see Figure 1. Thus in we use FEM to solve the equation

 ε∂ttE+∇×∇×E=0 in ΩFEM×(0,T),E(⋅,0)=f0(⋅), and ∂tE(⋅,0)=f1(⋅)% in ΩFEM,∂nE=−∂tE=g on ∂ΩFEM×(0,T),∇⋅(εE)=0 in ΩFEM. (15)

Here, is the restriction of the solution obtained by the FDM in to and therefore the test functions are not vanishing at the boundary and hence the term corresponding to in (9) will be appearing in the weak formulation.

To formulate the finite element method for (15) in , mimiking (9), we introduce the finite element space for each component of the electric field defined by

 WEh(ΩFEM):={w∈H1(ΩFEM):w|K∈P1(K),∀K∈Kh}.

Setting we define , , and to be the usual -interpolants of , , and , respectively, in . Then, similar to the FE scheme for , we get the following finite element scheme for : Given , , and , find such that

 (ε∂ttEh,v)+(∇Eh,∇v)+(∇⋅(εEh),∇⋅v)−(∇⋅Eh,∇⋅v)=⟨gh,v⟩∂ΩFEM,∀v∈WEh(ΩFEM)Eh(⋅,0)=f0h(⋅) and ∂tEh(⋅,0)=f1h(⋅) in ΩFEM. (16)

A corresponding fully discrete problem in reads as follows:

Given , , , , and ; find such that

 (Ek+1h,v)=(2Ekh,v)−(Ek−1h,v)−τ2(1/ε∇Ekh,∇v)−τ2(1/ε∇⋅(εEkh),∇⋅v)+τ2(1/ε∇⋅Ekh,∇⋅v)+τ2⟨gh/ε,v⟩∂ΩFEM,∀v∈WEh(ΩFEMT)E0h=f0hE1h=E0h+τf1h in ΩFEMT. (17)

Remark

Note that, in (15), Dirichlet boundary condition can be considered as well.

### 4.3 Fully discrete FE scheme for the electric field in ΩT

We expand the functions in terms of the standard continuous piecewise linear functions in space as

 (18)

where denote unknown coefficients at and the spatial mesh point . Then, substituting of (18) in (14), and setting , we obtain the linear system of equations:

 MEk+1h=2MEkh−MEk−1h−τ2G1Ekh−τ2G2Ekh+τ2G3Ekh+τ2M∂KEkh. (19)

Note that, unlike (14), now the contributions at boundary of the element appear in . Here, are the block mass matrices in space, are the block stiffness matrices in space, denote the nodal values of , and is a uniform time step. Now we define the mapping from the reference element onto such that and let be the piecewise linear local basis function on such that . Then, the explicit formulas for the entries in system of equations (19), for each element , can be written as:

 MKi,j=(φi(x)∘FK,φj(x)∘FK)K,M∂Ki,j=⟨1ε∂nφi(x)∘FK,φj(x)∘FK⟩∂K,G1Ki,j=(1ε∇φi∘FK,∇φj∘FK)K,G2Ki,j=(1ε∇⋅(εφi)∘FK,∇⋅φj∘FK)K,G3Ki,j=(1ε∇⋅φi∘FK,∇⋅φj∘FK)K, (20)

where , and , denote the , and , scalar products on and , respectively. Note that here is only the part of the boundary of element that lies at .

To obtain fully explicit scheme we approximate with the lumped mass matrix , (see [14, 18, 6] for the details corresponding to the Maxwell’s system (2)). Next, we multiply (20) by , and get the following explicit, fully discrete method in :

 Ek+1h=2Ekh−Ek−1h−τ2(ML)−1G1Ekh−τ2(ML)−1G2Ekh+τ2(ML)−1G3Ekh+τ2(ML)−1M∂KEkh. (21)

### 4.4 Fully discrete scheme for the electric field in ΩFEMT

As in the fully discrete FE scheme (19) in , we obtain fully discrete FE scheme in in the domain decomposition setting: Expanding the functions of via the continuous piecewise linear functions in space as in (18), and then substituting them in (17), (with , and ), we get the linear system of equations:

 MEk+1h=2MEkh−MEk−1h−τ2G1Ekh−τ2G2Ekh+τ2G3Ekh+τ2S∂K. (22)

Here, is the block mass matrices in space, restricted to , otherwise the same as in (19), are the block stiffness matrices in space as in (19),

is the assembled load vector ,

denote the nodal values of , is the time step. All quantities are for . Defining the mapping for the reference element in the mesh generated in as in the previous section, the formulas for entries of all matrices in the system (22) are the same as those in (20), and the entries of load vector are computed as

 SKj=(ghε,φj∘FK)∂K. (23)

Again, approximating with the lumped mass matrix , we obtain the following fully explicit scheme:

 Ek+1h=2Ekh−Ek−1h−τ2(ML)−1G1Ekh−τ2(ML)−1G2Ekh+τ2(ML)−1G3Ekh+τ2(ML)−1S∂K. (24)

### 4.5 Finite difference formulation

We recall now that from conditions (3) it follows that in the function . This means that in for the model problem (2) the forward problem will be

 ∂2E∂t2−ΔE=0 in  ΩFDM×(0,T), (25) E(x,0)=f0(x),   Et(x,0)=f1(x) in  ΩFDM, (26) E=0 on ∂Ω×(0,T), (27) ∂nE=∂nEFEM on ∂ΩFDM×(0,T). (28)

Using standard finite difference discretization of the equation (25) in we obtain the following explicit scheme for the solution of the forward problem:

 Ek+1l,j,m=τ2ΔEkl,j,m+2Ekl,j,m−Ek−1l,j,m. (29)

In the, system of, equations above, is the solution at the time iteration at the discrete point , is the time step, and is the discrete Laplacian.

Note that, in (29), the Dirichlet boundary consitions can be considered as well.

## 5 Stability

In this section we derive stability estimates for the semi-discrete approximations. For stability in these estimates are extensions of the stability approach derived for the continuous problem in [4]. As for the stability in we get slightly different norms involving contributions corresponding to the reflexive boundary: . We use discrete version of a triple norm induced by the weak variational formulation of (6), where we use the relation (10) (which is not necessary in the continuous case where , however, in general ):

Find such that

 (ε∂ttE,v)+(∇E,∇v)+((∇ε)⋅E,∇⋅v)+((ε−1)∇⋅E,∇⋅v)=0,∀v% ∈WE0(Ω)E(⋅,0)=f0(⋅) and ∂tE(⋅,0)=f1(⋅) in Ω. (30)

Remark

In general, in non-divergent free case, the bilinear form induced by (30) is not coercive. Further -conforming finite element may result in spurious solutions. A remedy is through modifying the equation by adding a gauge constrain of Coulomb-type, see, e.g. [21] and [23]. This is supplied by the ”zero”-term: , in (6), which we add in the continuous variational formulation in (9). This, however, is not necessarily true in the discrete forms, e.g. in (16), where most likely Taking in (30), (we used the boundary condition on ), yields

 (ε∂ttE,∂tE)+(∇E,∇∂tE)+((∇ε)⋅E,∇⋅∂tE)+((ε−1)∇⋅E,∇⋅∂tE)≡0, (31)

which, due to the fact that is independent of , can be rewritten as

 (32)

Proposition [4]

Let be a bounded domain with piecewise linear boundary . Then, the equation (6) has a unique solution . Further Let and , then there is a constant such that, , the following stability estimate holds true

 (33)
###### Proof.

The estimate (33) is proved in [4], Theorem 4.1 by setting and . Integrating (33) over the time interval we get the desired result. We omit the details. ∎

Below we translate this stability to the semi-discrete problem.

### 5.1 Stability estimate for the semi-discrete problem in Ω

The stability for the semi-discrete problem in is basically as in the continuous case above where all :s are replaced by with some relevant assumptions in the discrete data, viz.

Lemma

Assume that the interpolants of the data and : and satisfy the regularity conditions and , then for each

 |||Eh|||2ε(t)≤Ctε(||f1,h||2ε+||∇f0,h||2+||f0,h||2+||∇⋅f0,h||2ε−1). (34)

where

 |||Eh|||2ε(t):=||∂tEh||2ε(t)+||∇Eh||2(t)+||∇⋅Eh||2ε−1(t). (35)

### 5.2 Stability of the semi-discrete problem in ΩFEM

The stability of the semi-discrete problem in , relying on the variational formulation (16), and due to the appearance of the data function , is slightly different from (34). We rewrite (16), in view of (10), and with as: Given , , and , find such that

 (ε∂ttEh,∂tEh)+(∇Eh,∇∂tEh)+((∇ε)⋅Eh,∇⋅∂tEh)+((ε−1)∇⋅Eh,∇⋅∂tEh)=⟨gh,∂tEh⟩∂ΩFEM. (36)

To deal with the -term we rewrite (36) in its original form as (9) for and with :

 (ε∂ttEh,∂tEh)+(∇Eh,∇∂tEh)+(∇⋅(εEh),∇⋅∂tEh)−(∇⋅Eh,∇⋅∂tEh)=⟨gh,∂tEh⟩∂ΩFEM.