A Schwarz Method for the Magnetotelluric Approximation of Maxwell's equations

TL;DR

This paper introduces a Schwarz iterative method for solving the magnetotelluric approximation of Maxwell's equations, with proven convergence and numerical validation, aiding subsurface resource detection.

Contribution

It develops a classical Schwarz method tailored for complex-valued magnetotelluric equations and proves its convergence using a novel maximum modulus principle.

Findings

Convergence of the Schwarz method is established.

Numerical experiments confirm the theoretical results.

The approach effectively models electromagnetic wave propagation in Earth's subsurface.

Abstract

The magnetotelluric approximation of the Maxwell's equations is used to model the propagation of low frequency electro-magnetic waves in the Earth's subsurface, with the purpose of reconstructing the presence of mineral or oil deposits. We propose a classical Schwarz method for solving this magnetotelluric approximation of the Maxwell equations, and prove its convergence using maximum principle techniques. This is not trivial, since solutions are complex valued, and we need a new result that the magnetotelluric approximations satisfy a maximum modulus principle for our proof. We illustrate our analysis with numerical experiments.

Equations

$$\frac{\partial B}{\partial t}+\nabla \times E=0,-\frac{\partial D}{\partial t}+\nabla \times H=J,$$

$$\nabla \times E=-i\omega \mu H,\nabla \times H=\sigma E+J^e.$$

$$E=(E_x,0,E_z)\text{and}H=(0,H_y,0).$$

$$\frac{\partial E_x}{\partial z}-\frac{\partial E_z}{\partial x}=-i\omega \mu H_y,$$

$$-\frac{\partial H_y}{\partial z}=\sigma E_x+J_x^e\text{or}{E}_x=-\frac{1}{\sigma }\frac{\partial H_y}{\partial z}-\frac{J_x^e}{\sigma },$$

$$\frac{\partial H_y}{\partial x}=\sigma E_z+J_z^e\text{or}{E}_z=\frac{1}{\sigma }\frac{\partial H_y}{\partial x}-\frac{J_z^e}{\sigma }.$$

$$-\frac{\partial }{\partial z}\left\{\frac{1}{\sigma }\frac{\partial H_y}{\partial z}\right\}-\frac{\partial }{\partial x}\left\{\frac{1}{\sigma }\frac{\partial H_y}{\partial x}\right\}+i\omega \mu H_y=\frac{\partial }{\partial z}\left\{\frac{J_x^e}{\sigma }\right\}-\frac{\partial }{\partial x}\left\{\frac{J_z^e}{\sigma }\right\}.$$

$$\Delta u-i\omega u=f\quad\mbox{in $\Omega$}\;,u=0\quad\mbox{on $\partial\Omega$\;.}$$

$$\begin{array}[]{rcll}\Delta u^{n+1}_{j}-i\omega u^{n+1}_{j}&=&f&\mbox{in $\Omega_{j}$},\\ u^{n+1}_{j}&=&0&\mbox{on $\overline{\Omega}_{j}\cap\partial\Omega$},\\ u^{n+1}_{j}&=&u^{n}_{\text{glob}}&\mbox{on $\Gamma_{j}$},\end{array}$$

$$u_{\text{glob}}^{n+1}=\sum _{j=1}^J{\chi }_j{u}_j^{n+1}.$$

$$\Delta u-i\omega u=0\text{mbox}and\Delta \overline{u}+i\omega \overline{u}=0.$$

$$\Delta |u{|}^2=\Delta (u\overline{u})=\nabla (\overline{u}\nabla u+u\nabla \overline{u})=\nabla \overline{u}\nabla u+\overline{u}\Delta u+\nabla u\nabla \overline{u}+u\Delta \overline{u}=$$

$$=2|\nabla u{|}^2+i\omega |u{|}^2-i\omega |u{|}^2=2|\nabla u{|}^2\ge 0.$$

$$\underset{j=1,\dots ,J}{\max }||u-u_j^n|{|}_{{\Omega }_j}\le {\gamma }^n\underset{j=1,\dots ,J}{\max }||u-u_j^0|{|}_{{\Omega }_j},$$

$$||e_j^{n+1}|{|}_{{\Omega }_j}=||e_j^{n+1}|{|}_{{\Gamma }_j}=||\sum _{j^{\prime }=1}^J{\chi }_{j^{\prime }}{e}_{j^{\prime }}^n|{|}_{{\Gamma }_j}$$

$$\le \underset{j^{\prime }=1,\dots ,J}{\max }||e_{j^{\prime }}^n|{|}_{K_{j^{\prime }}}\le \gamma \underset{j^{\prime }=1,\dots ,J}{\max }||e_{j^{\prime }}^n|{|}_{{\Gamma }_{j^{\prime }}}=\gamma \underset{j^{\prime }=1,\dots ,J}{\max }||e_{j^{\prime }}^n|{|}_{{\Omega }_{j^{\prime }}},$$

$$\underset{j=1,\dots ,J}{\max }||e_j^{n+1}|{|}_{{\Omega }_j}\le \gamma \underset{j^{\prime }=1,\dots ,J}{\max }||e_{j^{\prime }}^n|{|}_{{\Omega }_{j^{\prime }}},$$

Full text

1 Introduction

Maxwell’s equations can be used to model the propagation of electro-magnetic waves in the subsurface of the Earth. The interaction of such waves with the material in the subsurface produces response waves, which carry information about the physical properties of the Earth’s subsurface, and their measurement allows geophysicists to detect the presence of mineral or oil deposits. Since such deposits are often found to be invariant with respect to one direction parallel to the Earth’s surface, the model can be reduced to a two dimensional complex partial differential equation. Following [19], the magnetotelluric approximation is derived from the full 3D Maxwell’s equations,

[TABLE]

in the quasi-static (i.e. long wavelength, low frequency) regime, which implies that $\frac{\partial D}{\partial t}$ in (1) is neglected. Assuming a time dependence of the form $e^{i\omega t}$, where $\omega$ is the pulsation of the wave, using Ohm’s law, $J=\sigma E+J^{e}$, where $J^{e}$ denotes some exterior current source, and the constitutive relation $B=\mu H$ where $\mu$ is assumed to be the permeability of free space, we obtain

[TABLE]

Assuming the plane-wave source of magnetotellurics, and a two-dimensional Earth structure such that $\sigma=\sigma(x,z)$, the electric and magnetic fields can be decomposed into two independent modes. For the TM-, or H-polarization, mode, we have

[TABLE]

Hence the first vector valued equation in (2) becomes a scalar equation,

[TABLE]

and the second vector valued equation in (2) gives two scalar equations,

[TABLE]

and

[TABLE]

Substituting (4) and (5) into (3) thus leads to a scalar equation for $H_{y}$,

[TABLE]

In geophysical applications the coefficient of conductivity $\sigma$ is in general a non-constant, piece-wise continuous function. We will assume however for simplicity of our first study here that $\sigma\equiv 1$. If we then set $u:=H_{y}$, assume homogeneous Dirichlet boundary conditions and denote by $f:=-\frac{\partial}{\partial z}\left\{\frac{J_{x}^{e}}{\sigma}\right\}+\frac{\partial}{\partial x}\left\{\frac{J_{z}^{e}}{\sigma}\right\}$, we obtain as the magnetotelluric approximation of the Maxwell equations (cf. equation (2.86) in [19])

[TABLE]

We further assume for simplicity that $\Omega$ is a domain with smooth boundary and that $f\in C^{\infty}(\Omega)\cap C(\overline{\Omega})$. The pulsation $\omega$ is assumed to be real and non-zero. Note that the solution $u$ of equation (7) could also represent a component of the electric field if the model had been derived in an analogous fashion from the TE mode.

We are interested in solving the magnetotelluric approximation (7) using Schwarz methods. The alternating Schwarz method, introduced by H.A. Schwarz in 1869 [17] to prove existence and uniqueness of solutions to Laplace’s equation on irregular domains, is the foundational idea of the field of domain decomposition, and has inspired work in both theoretical aspects and applications to all fields of science and engineering, see [9, 4] and references therein for more information about the historical context. Lions [15, 16] reconsidered the problem of the convergence of the method for the Poisson equation on more general configurations of overlapping subdomains. In his second paper [16], he followed the idea of Schwarz and proved convergence of the alternating Schwarz method using the maximum principle for harmonic functions. He also introduced a parallel variant of the Schwarz method, where all subdomain problems are solved simultaneously. Schwarz methods have also been introduced and studied for the original Maxwell equations (1), see [2, 1, 6, 5, 7, 8] and references therein; maximum principle arguments can not be used here for studying convergence. We show here that for the magnetotelluric approximation of Maxwell equations in (7), which also has complex solutions like the original Maxwell equations, the convergence of the parallel Schwarz method can be proved using a maximum modulus principle satisfied by complex solutions of (7).

2 Well-Posedness, Schwarz Method and Convergence

We start by establishing the well-posedness of the magnetotelluric approximation of the Maxwell equations in (7).

Theorem 1**.**

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded domain with smooth boundary. Assume that $f\in L^{2}(\Omega)$ and $\omega$ is a non-zero constant. Then the boundary value problem (7) has a unique solution $u\in H^{1}_{0}(\Omega)$, depending continuously on $f$.

Proof.

This result follows from a standard application of the Riesz Representation Theorem and the Lax-Milgram Lemma. ∎

We now decompose the domain $\Omega\subset\mathbb{R}^{2}$ first into non-overlapping subdomains, and then enlarge each subdomain by a layer of positive width to obtain the overlapping subdomains $\Omega_{j}$, for $j=1,...J$, leading to a strongly overlapping subdomain decomposition of $\Omega$. An example is shown in Figure 1, where the non-overlapping decomposition is indicated by the dashed lines, see also [10].

For such strongly overlapping decompositions, one can define a smooth partition of unity $\{\chi_{j}\}_{j=1}^{J}$ subordinated to the open covering $\{\Omega_{j}\}_{j=1}^{J}$, such that the support of $\chi_{j}$ is a set $K_{j}$ contained in the open subdomain $\Omega_{j}$ for each $j=1,2,\ldots,J$, see [18, Theorem 15, Chapter 2]. The assumption of a strongly overlapping decomposition is not strictly necessary to use maximum principle arguments, see for example [11, 12], which contain even accurate convergence estimates, but we make it here since it simplifies the new and first application of the maximum modulus principle (via Corollary 4) for studying Schwarz methods for equations with complex valued solutions. For each $\Omega_{j}$, we denote by $\Gamma_{j}$ the portion of $\partial\Omega_{j}$ in the interior of $\Omega$.

The parallel Schwarz method for such a multi-subdomain decomposition starts with a global initial guess for the solution of (7), $u_{\text{glob}}^{0}\in C^{2}(\Omega)\cap C^{0}(\overline{\Omega})$ (less regularity would also be possible, because of the regularization provided by the equation). If at the step $n$ of the parallel Schwarz method the global approximation $u_{\text{glob}}^{n}$ has been constructed, and $u_{\text{glob}}^{n}\in C^{2}(\Omega)\cap C^{0}(\overline{\Omega})$, then it constructs the next global approximation at step $n+1$ by solving, for $j=1,\cdots,J$, the Dirichlet problems

[TABLE]

and then defining the $(n+1)^{th}$ global iterate by using the partition of unity,

[TABLE]

Since the initial guess $u_{\text{glob}}^{0}$ is smooth, by induction it follows that $u_{\text{glob}}^{n+1}\in C^{2}(\Omega)\cap C^{0}(\overline{\Omega})$. This fact allows us to use the classical (i.e. non-variational ) formulation of the maximum modulus principle.

Definition 2**.**

A real valued function $v$ of class $C^{2}(\Omega)$ is said to be subharmonic if $\Delta v(x)\geq 0$, $\forall x\in\Omega$, and strictly subharmonic if $\Delta v(x)>0$, $\forall x\in\Omega$.

Note that the above definition is not the most general one, but it is suitable for the purposes of our paper. The property that we will use to prove the convergence of the parallel Schwarz method is the well-known maximum principle, which is the content of the next theorem (see [13], Theorem J-7).

Theorem 3**.**

Let $v\in C^{2}(\Omega)\cap C(\overline{\Omega})$ be a non-constant subharmonic function. Let $O\subset\Omega$ be a proper open subset. Then $v$ satisfies the strong maximum principle, namely $\max_{O}v<\max_{\partial\Omega}v$.

By continuity one can immediately deduce the following corollary, which contains the key estimate for proving the convergence of the parallel Schwarz method.

Corollary 4**.**

Let $u\in C^{2}(\Omega)\cap C^{0}(\overline{\Omega})$ be a non-constant subharmonic function, and $K$ be a closed subset of $\Omega$. Then there exists a constant $\gamma\in[0,1)$ such that $\max_{K}u<\gamma\max_{\partial\Omega}u$.

Proof.

It follows by observing that since $K$ is closed, it has a positive distance from the boundary $\partial\Omega$ of the open domain $\Omega$. ∎

Since the solution of the magnetotelluric approximation (7) of Maxwell’s equation has complex valued solutions, it is not directly possible to use the maximum principle result in Corollary 4 for proving convergence of the associated Schwarz method (8)-(9). The key additional ingredient is to prove the following property on the modulus of solutions of the magnetotelluric approximation:

Theorem 5**.**

Let $u\in C^{2}(\Omega)\cap C(\overline{\Omega})$ be a non-zero solution of the homogeneous form of equation (7). Then $|u|^{2}$ is a non-constant subharmonic function.

Proof.

Taking the complex conjugate of the partial differential equation (7) with $f=0$, gives a pair of equations,

[TABLE]

Hence we can compute

[TABLE]

Therefore we have shown that $|u|^{2}$ is subharmonic. If $|u|^{2}$ is constant, the same calculations show that $\nabla u\equiv 0$, which implies that $u$ is a constant solution, hence it must be identically equal to zero, since the equation $\Delta u-i\omega u=0$ has no constant non-zero solutions. ∎

We can now prove our main result, namely the convergence of the parallel Schwarz method for solving the magnetotelluric approximation of Maxwell’s equation in the infinity norm, which we denote by $||\cdot||_{S}$ for any function on a set $S$.

Theorem 6**.**

The parallel Schwarz method (8)-(9) for the magnetotelluric approximation (7) of Maxwell’s equations is convergent and satisfies the error estimate

[TABLE]

where $u$ denotes the global solution of problem (7) and $u_{j}^{n}$ the approximations from the parallel Schwarz method (8)-(9), and the constant $\gamma<1$ comes from Corollary 4.

Proof.

For $j=1\cdots J$, let $K_{j}\subset\Omega_{j}$ be the support of the partition of unity function $\chi_{j}$, and let $e_{j}^{n}:=u-u_{j}^{n}$ be the error. Then $e_{j}^{n}$ is solution of the homogeneous equation $\Delta e_{j}^{n}-i\omega e_{j}^{n}=0$, and hence by Theorem 5 its modulus is a subharmonic function, and thus by Theorem 3, the modulus of the error error $|e^{n}_{j}|$ satisfies the strong maximum principle. We can then estimate on each subdomain $\Omega_{j}$

[TABLE]

where $\gamma\in[0,1)$ is the maximum of the factor introduced in Corollary 4 over all $\Omega_{j}$ and corresponding $K_{j}$. Since this holds for all $j$, we can take the maximum on the left and obtain

[TABLE]

which proves by induction (10). ∎

Remark 7**.**

The convergence factor $\gamma<1$ is not quantified in Theorem 6, since Corollary 4 does not provide a method to estimate the constant $\gamma$ in the generality of the decomposition we used, but such an estimate is possible for specific decompositions, see for example [11, 12].

Remark 8**.**

In [16], Lions proved the convergence of the classical Schwarz method for the Poisson equation with Dirichlet boundary conditions using a method that does not use the maximum principle. His remarkable proof is based on the method of orthogonal projections, and relies on the fact that the bilinear form associated with the weak formulation of the Poisson equation is an inner product in the solution space $H^{1}_{0}(\Omega)$. We do not see how this method can be extended to prove convergence of the classical Schwarz method applied to the magnetotelluric approximation of the Maxwell equations. In our case in fact the bilinear form associated to the weak formulation of (7) of the global problem is not an inner product, as it fails to be symmetric and positive-definite.

3 Numerical examples

We now present two numerical experiments. The simulations are computed on a domain $\Omega$ that consists of two squares $\Omega_{1}$ and $\Omega_{2}$, each of unit size $1\times 1$. The discretization for each square consists of a uniform grid of $30\times 30$ points. The overlap is along a vertical strip whose width is specified by the number of grid points, denoted by $d$.

We first compute the error $e_{j}^{n}:=u-u_{j}^{n}$, as used in the proof of Theorem 6. In Figure 2

we show, from left to right, the modulus of the error on the left subdomain for iteration $n=1$, $n=5$ and $n=15$, for an overlap of $d=6$ horizontal grid points. We chose $\omega=1$, and the initial error was produced by generating random values uniformly distributed on the range $[0,1]$. Note how the modulus of the error clearly satisfies the maximum principle. In Figure 3,

we plot the dependence of the interface residual in the 2-norm on the iteration number, for three different overlap sizes $d=2,4,6$. As expected, the performance of the algorithm improves as we increase the size of the overlap, since increasing the overlap improves the coefficient $\gamma$ in Corollary 4 which is the key quantity governing the convergence of the parallel Schwarz method.

4 Conclusion

We showed in this paper that even though the solutions of the magnetotelluric approximation of Maxwell’s equations are complex valued, maximum principle arguments can be used to prove convergence of a parallel Schwarz method. The main new ingredient is a maximum modulus principle which is satisfied by the solutions of the magnetotelluric approximation. In a forthcoming paper, we will analyze the convergence rate of the parallel Schwarz method via Fourier analysis, and we will also introduce more efficient transmission conditions of Robin (or higher-order) type at the interfaces between the subdomains, which leads to optimized Schwarz methods, see [14, 3] and references therein.

Acknowledgement We would like to thank Dr. Hormoz Jandahari and Dr. Colin Farquharson for the insightful explanation of the physics beyond the model discussed in this manuscript. The research of the first author was funded, in part, by the Canada Research Chairs program, the NSERC Discovery Grant program and the IgniteR&D program of the Research and Development Corporation of Newfoundland and Labrador (RDC).