Domain Decomposition with Neural Network Interface Approximations for
  time-harmonic Maxwell's equations with different wave numbers
In this work, we consider the time-harmonic Maxwell's equations and their
numerical solution with a domain decomposition method. As an innovative
feature, we propose a feedforward neural network-enhanced approximation of the
interface conditions between the subdomains. The advantage is that the
interface condition can be updated without recomputing the Maxwell system at
each step. The main part consists of a detailed description of the construction
of the neural network for domain decomposition and the training process. To
substantiate this proof of concept, we investigate a few subdomains in some
numerical experiments with low frequencies. Therein the new approach is
compared to a classical domain decomposition method. Moreover, we highlight
current challenges of training and testing with different wave numbers and we
provide information on the behaviour of the neural-network, such as convergence
of the loss function, and different activation functions.
T. Knoke, S. Kinnewig, S. Beuchler, A. Demircan, U. Morgner, T. Wick
Domain Decomposition with Neural Network Interface Approximations for time-harmonic Maxwell's equations with different wave numbers

###### Abstract

In this work, we consider the time-harmonic Maxwell's equations and their numerical solution with a domain decomposition method. As an innovative feature, we propose a feedforward neural network-enhanced approximation of the interface conditions between the subdomains. The advantage is that the interface condition can be updated without recomputing the Maxwell system at each step. The main part consists of a detailed description of the construction of the neural network for domain decomposition and the training process. To substantiate this proof of concept, we investigate a few subdomains in some numerical experiments with low frequencies. Therein the new approach is compared to a classical domain decomposition method. Moreover, we highlight current challenges of training and testing with different wave numbers and we provide information on the behaviour of the neural-network, such as convergence of the loss function, and different activation functions.

**Keywords:** Time-Harmonic Maxwell's Equations, Machine Learning, Feedforward Neural Network, Domain Decomposition Method.

## 1 Introduction

The Maxwell's equations for describing electro-magnetic phenomena are of great interest in current research fields, such as optics. One present example of employing Maxwell's equations can be found in the Cluster of Excellence PhoenixD (Photonics Optics Engineering Innovation Across Disciplines)1 at the Leibniz University Hannover, in which modern methods for optics simulations are being developed. Therein, one focus is on the efficient and accurate calculation of light distribution in an optical material to design optical devices on the micro- and nanoscale [35, 26]. In comparison to other partial differential equations, such as in solid mechanics or fluid flow, the Maxwell's equations have some peculiarities such as the curl operator, which has in two-dimensional problems, a one-dimensional image, but in three-dimensional problems, it has a three-dimensional image. Moreover, the requirements for the discretization and definiteness of the final linear system are specific. In more detail, in numerical mathematics, Maxwell's equations are of interest because of their specific mathematical structures [27, 9, 25, 33], requirements for finite elements [27, 25, 28, 31, 21, 7, 29], their numerical solution [18, 16, 14, 10] as well as postprocessing such as a posteriori error control and adaptivity [34, 6]. As their numerical solution is challenging due to their ill-posed nature, e.g., [4], one must apply suitable techniques. The most prominent approach in the literature is based on domain decomposition (DD) techniques [36, 10]. The geometric multigrid solver developed by Hiptmair [18] can only be applied to the problem in the time domain (i.e., the well-posed problem).
In this work, we concentrate on the numerical solution using a domain decomposition method. Specifically, our starting point is the method developed in [4], based on ideas from [11], and which was realized in the modern open-source finite element library deal.II [1, 2]. The domain decomposition method's crucial point is the interface operator derivation [11]. Our main objective in the current work is to design a proof of concept to approximate the interface operator with the help of a feedforward neural network (NN) [5, 17, 20]. We carefully derive the governing algorithms and focus on a two-domain problem to study our new approach's mechanism and performance. Implementation-wise, the previously mentioned deal.II library (in C++) is coupled to the PyTorch (in python) [32] library, which is one of the standard packages for neural network computations. Our main aim is to showcase that our approach is feasible and can be a point of departure for further future extensions. We notice that the current work is an extension of the conference proceedings paper [22] with more mathematical and algorithmic details, and different numerical tests, specifically the studies on different wave numbers and comparison of two NN activation functions.

The outline of this work is as follows: In Section 2, we introduce the time-harmonic Maxwell's equations and our notation. Next, in Section 3, domain decomposition and neural network approximations are introduced. Afterwards, we address in detail the training process in Section 4. In Section 5, some numerical tests demonstrate our proof of concept. Our work is summarized in Section 6.

## 2 Equations

For the sake of simplicity, we only consider the two-dimensional time-harmonic Maxwell's equations. In the following, we will introduce these equations in detail.

### Fundamental operators

To comprehensively describe the problem, we introduce the basic operators needed to describe two-dimensional electro-magnetic problems. Therefore, let us assume a scalar function \(\phi:\mathbb{R}\rightarrow\mathbb{R}\) and \(\vec{v}\in\mathbb{R}^{2}\) to be a two-dimensional vector. Then the gradient of \(\phi\) is given by \(\nabla\phi=\left(\frac{\partial\phi}{\partial x_{1}},\ \frac{\partial\phi}{ \partial x_{2}}\right),\) and the divergence of \(v\) is given by \(\mathrm{div}(v)\coloneqq\nabla\cdot v\coloneqq\sum_{i=1}^{2}\frac{\partial v_ {i}}{\partial x_{i}}.\) Next, \(a\cdot b=(a_{1},a_{2})^{T}\cdot(b_{1},b_{2})^{T}=a_{1}b_{1}+a_{2}b_{2}\) denotes the scalar product. We can furthermore write down the description of the two-dimensional curl operator

\[\mathrm{curl}(\vec{v})=\frac{\partial v_{2}}{\partial x_{1}}-\frac{\partial v _{1}}{\partial x_{2}}, \tag{1}\]

and the curl operator applied to a scalar function

\[\underline{\mathrm{curl}}(\phi)=\left(\begin{array}{c}\frac{\partial\phi} {\partial x_{2}}\\ -\frac{\partial\phi}{\partial x_{1}}\end{array}\right). \tag{2}\]

### Time-harmonic Maxwell's equations

Let \(\Omega\subset\mathbb{R}^{2}\) be a bounded domain with sufficiently smooth boundary \(\Gamma\). The latter is partitioned into \(\Gamma=\Gamma^{\infty}\cup\Gamma^{\mathrm{inc}}\). The main governing function space is defined as

\[H(\mathrm{curl},\Omega):=\{v\in\mathcal{L}^{2}(\Omega)\ |\ \mathrm{curl}(v)\in \mathcal{L}^{2}(\Omega)\},\]

where \(\mathcal{L}^{2}(\Omega)\) is the well-known space of square-integrable functions in the Lebesgue sense. In order to define boundary conditions, we introduce the traces

\[\gamma^{t}:H(\mathrm{curl},\Omega)\to H_{\times}^{-1/2}( \mathrm{div},\Gamma),\] \[\gamma^{T}:H(\mathrm{curl},\Omega)\to H_{\times}^{-1/2}( \mathrm{curl},\Gamma),\]

which are defined by

\[\gamma^{t}\left(\phi\right)=\left(\begin{array}{c}\phi\ n_{2}\\ -\phi\ n_{1}\end{array}\right)\quad\text{and}\quad\gamma^{T}\left(v\right)=v- \left(n\cdot v\right)\cdot n,\]
where \(n\in\mathbb{R}^{2}\) is the normal vector of \(\Omega\), \(H_{\times}^{-1/2}(\operatorname{div},\Gamma):=\{v\in H^{-1/2}(\Gamma)\mid v\cdot n =0,\;\operatorname{div}_{\Gamma}v\in H^{-1/2}(\Gamma)\}\) is the space of well-defined surface divergence fields and \(H(\operatorname{curl},\Gamma):=\{v\in H^{-1/2}(\Gamma)\mid v\cdot n=0,\; \operatorname{curl}_{\Gamma}\left(v\right)\in H^{-1/2}(\Gamma)\}\) is the space of well-defined surface curls, see [27, Chapter 3.4]. In the following, we first state the strong form of the system. The time-harmonic Maxwell's equations are then defined as follows: Find the electric field \(E:\Omega\to\mathbb{C}^{2}\) such that

\[\left\{\begin{array}{ll}\underline{\operatorname{curl}}\left(\mu^{-1} \operatorname{curl}\left(\vec{E}\right)\right)-\varepsilon\omega^{2}\vec{E}& =\vec{0}\qquad\text{ in }\Omega\\ \mu^{-1}\gamma^{t}\left(\operatorname{curl}\left(\vec{E}\right)\right)-i \kappa\omega\gamma^{T}\left(\vec{E}\right)&=\vec{0}\qquad\text{ on }\Gamma^{\infty}\\ \gamma^{T}\left(\vec{E}\right)&=\vec{E}^{\text{inc}}\quad\text{ on }\Gamma^{\text{inc}},\end{array}\right. \tag{3}\]

where \(\vec{E}^{\text{inc}}:\mathbb{R}^{2}\to\mathbb{C}^{2}\) is some given incident electric field, \(\mu\in\mathbb{R}^{+}\) is the relative magnetic permeability, \(\kappa=\sqrt{\varepsilon}\), \(\varepsilon\in\mathbb{C}\) relative permittivity, \(\omega=\frac{2\pi}{\lambda}\) is the wave number and \(\lambda\in\mathbb{R}^{+}\) is the wave length and \(i\) denotes the imaginary number. System (3), as well as its weak form, is called time-harmonic, because the time dependence can be expressed by \(e^{i\omega\tau}\), where \(\tau\geq 0\) denotes the time.

### Weak formulation

In this subsection, we derive the weak form. This is the starting point for a finite element method (FEM) discretization. For the derivation, we first begin by rewriting the curl product with the help of integration by parts:

\[\int_{\Omega}\underline{\operatorname{curl}}\left(\phi\right)\cdot\vec{u}\; \mathsf{d}x=\int_{\Omega}\phi\operatorname{curl}\left(u\right)\;\mathsf{d}x+ \int_{\partial\Omega}\gamma^{t}(\phi)\cdot u\;\mathsf{d}s, \tag{4}\]

see for instance [15, 27]. We want to derive the weak formulation from the strong formulation (3) in the following:

\[\int_{\Omega}\underline{\operatorname{curl}}\left(\mu^{-1} \operatorname{curl}\left(\vec{E}\right)\right)\cdot\vec{\varphi}\;\mathsf{d}x- \varepsilon\omega^{2}\int_{\Omega}\vec{E}\cdot\vec{\varphi}\;\mathsf{d}x= \vec{0},\] \[\stackrel{{\eqref{eq:weakform}}}{{\Rightarrow}}\int_{ \Omega}\mu^{-1}\operatorname{curl}\left(\vec{E}\right)\operatorname{curl} \left(\vec{\varphi}\right)\;\mathsf{d}x-\varepsilon\omega^{2}\int_{\Omega} \vec{E}\cdot\vec{\varphi}\;\mathsf{d}x+\int_{\partial\Omega}\mu^{-1}\gamma^{t} \left(\operatorname{curl}\left(\vec{E}\right)\right)\cdot\vec{\varphi}\; \mathsf{d}s=\vec{0}. \tag{5}\]

By applying the definition of the boundaries \(\Gamma^{\infty}\) and \(\Gamma^{\text{inc}}\) from equation (3) to equation (5), we obtain the weak formulation of the time-harmonic Maxwell's equations. Find \(\vec{E}\in H(curl,\Omega)\) such that for all \(\vec{\varphi}\in H(curl,\Omega)\)

\[\int_{\Omega}\left(\mu^{-1}\operatorname{curl}\left(\vec{E} \right)\operatorname{curl}\left(\vec{\varphi}\right)-\varepsilon\omega^{2}\vec {E}\cdot\vec{\varphi}\right)\;\mathsf{d}x+ i\kappa\omega\int_{\Gamma^{\infty}}\gamma^{T}\left( \vec{E}\right)\cdot\gamma^{T}\left(\vec{\varphi}\right)\;\mathsf{d}s\] \[= \int_{\Gamma^{\text{inc}}}\gamma^{T}\left(\vec{E}^{\text{inc}} \right)\cdot\gamma^{T}\left(\vec{\varphi}\right)\;\mathsf{d}s. \tag{6}\]

### Two-dimensional Nedelec elements

For the implementation with the help of a Galerkin finite element method (FEM), we need the discrete weak form. Based on the De-Rham cohomology, we must choose our basis functions out of the Nedelec space \(V_{h}\). Therefore, we want to introduce the definition of the space \(V_{h}\) in the following, based on the formalism introduced by Zaglmayr [37, Chapter 5.2].

As a suitable polynomial basis, we introduce the integrated Legendre polynomials. Let \(x\in[-1,1]\). The following recursive formula defines the integrated Legendre polynomials:

\[L_{1}(x) =x,\] \[L_{2}(x) =\tfrac{1}{2}\left(x^{2}-1\right), \tag{7}\] \[(n+1)L_{n+1}(x) =(2n-1)xL_{n}(x)-(n-2)L_{n-1}(x),\quad\text{ for }n\geq 2.\]

Let us choose the quadrilateral reference element as \(Q=[0,1]\times[0,1]\).
We continue by defining the set of all edges with local edge-ordering, where, see Figure 1. We denote the cell itself with local vertex-ordering. The polynomial order is given by.

We continue by defining the set of all edges with local edge-ordering, where. We denote the cell itself with local vertex-ordering. The polynomial order is given by.

With the help of these basis functions, we define the two-dimensional Nedelec space

(8)

where is the space of the _lowest-order Nedelec_ function, is the space of the _edge-bubbles_ and is the space of the _cell-bubbles_. All basis functions on one element with barycentric coordinates are displayed in Figure 2. Visualizations of some basis functions are displayed in Figure 3. The description of is still not complete, so far we only described, with as previously defined. It remains to introduce the Piola transformation, which is used to transform the reference element to any given physical element, see Monk [27] (Lemma 3.57, Corollary 3.58).

### Discrete weak formulation

We have gathered everything to write down the discrete weak formulation of the time-harmonic Maxwell's equations. We obtain the discrete weak formulation by applying the Galerkin method

Figure 1: Left: Vertex and edge ordering on the reference cell, right: parametrisation of the reference cell.
to the equation (6). Find \(E_{h}\in V_{h}(\Omega)\) such that

\[\int_{\Omega}\left(\mu^{-1}\operatorname{curl}\left(\vec{E}_{h} \right)\operatorname{curl}\left(\vec{\varphi}_{h}\right)-\varepsilon\omega^{2} \vec{E}_{h}\cdot\vec{\varphi}_{h}\right)\ \mathsf{d}x+ i\kappa\omega\int_{\Gamma^{\infty}}\gamma^{T}\left(\vec{E}_{h} \right)\cdot\gamma^{T}\left(\vec{\varphi}_{h}\right)\ \mathsf{d}s\] \[=\int_{\Gamma^{\mathrm{inc}}}\gamma^{T}\left(\vec{E}^{\mathrm{ inc}}\right)\cdot\gamma^{T}\left(\vec{\varphi}_{h}\right)\ \mathsf{d}s\ \ \forall\varphi_{h}\in V_{h}(\Omega). \tag{9}\]

## 3 Numerical approach

In this section, we first describe domain decomposition and afterwards the neural network approximation. In the latter, we also outline how to replace the interface operator by the neural network.

### Domain decomposition

Since the solution of Maxwell's equation system (3) is challenging, as already outlined in the introduction, we apply a non-overlapping domain composition method (DDM)[36] in which the domain is divided into subdomains as follows

\[\overline{\Omega}=\bigcup_{i=0}^{n_{\mathrm{dom}}}\overline{ \Omega}_{i}\quad\text{with}\] \[\Omega_{i}\cap\Omega_{j}=\varnothing\quad\forall i\neq j,\]

where \(n_{\mathrm{dom}}+1\) is the number of subdomains. In such a way, every subdomain \(\Omega_{i}\) becomes small enough so that we can handle it with a direct solver. The global solution of the electric field \(E\) is computed via an iterative method, where we solve the time-harmonic Maxwell's equations on each subdomain with suitable interface conditions between the different subdomains. Thus, we obtain a solution \(E_{i}^{k}\) for

Figure 3: Plots of basis functions on \((0,1)^{2}\): low order edge function (left above), high order edge based basis function for \(p=2\) (right above) to the edge \(\mathcal{E}_{0}\), high order cell based basis functions for \(p=2\) of type 1 and 2 (below).
every subdomain \(\Omega_{i}\), where \(k\) denotes the \(k\)-th iteration step. The initial interface condition is given by

\[g_{ji}^{k=0}:=-\mu^{-1}\gamma_{i}^{t}\left(\operatorname{curl}\left(E_{i}^{k=0} \right)\right)-i\kappa S\left(\gamma_{i}^{T}\left(E_{i}^{k=0}\right)\right)=0, \tag{10}\]

where \(S\) describes the interface operator, \(i\) is the index of the current domain, and \(j\) is the index of the neighbouring domain [11]. Afterwards, the electric-field \(E_{i}^{k+1}\) is computed at each step by solving the following system

\[\left\{\begin{array}{ll}\operatorname{curl}\left(\mu^{-1} \operatorname{curl}\left(E_{i}^{k+1}\right)\right)-\omega^{2}\varepsilon E_{i }^{k+1}&=0&\text{in }\Omega_{i},\\ \mu^{-1}\gamma_{i}^{t}\left(\operatorname{curl}\left(E_{i}^{k+1}\right) \right)-i\omega\kappa\gamma_{i}^{T}\left(E_{i}^{k+1}\right)&=0&\text{on }\Gamma_{i}^{ \infty},\\ \gamma_{i}^{T}\left(E_{i}^{k+1}\right)&=\gamma_{i}^{T}\left(E_{i}^{\text{inc} }\right)&\text{on }\Gamma_{i}^{\text{inc}},\\ \mu^{-1}S\left(\gamma_{i}^{t}\left(\operatorname{curl}\left(E_{i}^{k+1} \right)\right)\right)-i\omega\kappa\gamma_{i}^{T}\left(E_{i}^{k+1}\right)&=g_{ ji}^{k}&\text{on }\Sigma_{ij},\end{array}\right. \tag{11}\]

where \(\Sigma_{ij}=\Sigma_{ji}:=\partial\Omega_{i}\cap\partial\Omega_{j}\) denotes the interface of two neighbouring elements and the interface condition is updated by

\[g_{ji}^{k+1}=-\mu^{-1}\gamma_{i}^{t}\left(\operatorname{curl}\left(E_{i}^{k+1 }\right)\right)-i\kappa S\left(\gamma_{i}^{T}\left(E_{i}^{k+1}\right)\right)=- g_{ij}^{k}-2i\kappa S\left(\gamma_{i}^{T}\left(E_{i}^{k+1}\right)\right). \tag{12}\]

In case of success we obtain \(\lim_{k\to\infty}E_{i}^{k}=E|_{\Omega_{i}}\), but this convergence depends strongly on the chosen interface operator \(S\) (see [10, 11]). The implementation of this approach into deal.II was done in [4].

### Our new approach: Neural network approximation of \(S\)

Since the computation of a good approximation of \(S\) is challenging, we examine a new approach in which we attempt to approximate this operator with the help of a neural network (NN). For a first proof of concept, we choose a prototype example and explore whether an NN can approximate the interface values. As it is not feasible to compute the exact interface operator \(S\), we aim to compute \(g_{ij}^{k+l},\ l>0\) with an NN, using \(g_{ij}^{k}\) and \(E_{i}^{k+1}\) as input. Another benefit of this approach is that we can quickly generate a training data set from a classical domain decomposition method, as described in Section 4.4. We choose \(S=\mathds{1}\) for simplicity inside our classical domain decomposition method. Hence, the advantage of this approach is that one can update the interface condition without recomputing the system (11) at each step, raising the hope of reducing the computational cost.

## 4 Neural network training

The first step in neural network approximations is the training process, which is described in this section. Besides the mathematical realization, we also need to choose the software libraries. We utilize deal.II [2] to discretize the time-harmonic Maxwell's equations with the finite element method. The neural network is trained with PyTorch [32]. The exchange of information between the results of the deal.II code and the PyTorch code take place via the hard disk.

### Basic definitions

First of all, we give a short definition of the neural network type employed in this work, and we introduce the basic parameters. Further information can be found in [5, 8, 24, 3, 30]. The following notation and descriptions of this subsection are mainly based on [23].

**Definition 4.1** (Artificial neuron).: _An (artificial) neuron (also known as unit [12], [5][Section 5.1]) \(u\) is a tuple of the form \((\mathfrak{x},\mathfrak{w},\sigma)\). The components have the following meanings:_

* \(\mathfrak{x}=(\mathfrak{x}_{0},\ldots,\mathfrak{x}_{n})\in\mathbb{R}^{n+1}\) _is the input vector. It contains the information, that the neuron receives._
* \(\mathfrak{w}=(\mathfrak{w}_{0},\ldots,\mathfrak{u}_{n})\in\mathbb{R}^{n+1}\) _is the weight vector, which determines the influence of the individual input information on the output of the neuron. Later,_ \(\mathfrak{w}\) _denotes the weight vector of all neurons._
* \(\sigma:\mathbb{R}\to\mathbb{R}\)_, with_ \(\sigma=\sum_{i=0}^{n}\mathbbm{x}_{i}\mathbbm{u}_{i}\mapsto a\) _is the activation function. It determines the so-called activation level_ \(a\) _from the input and the weights, which represent the output of the neuron._

**Definition 4.2** (Neural network).: _An (artificial) neural (feedforward-) network is a set of neurons \(U\) with a disjoint decomposition \(U=U_{0}\dot{\cup}\ldots\dot{\cup}U_{l}\). The partition sets \(U_{\mathsf{k}},\mathsf{k}=0,\ldots,l\) are called layers. Here, \(U_{0}\) is the input layer. It contains the neurons that receive information from outside. Moreover, \(U_{l}\) is the output layer with the neurons that return the output. Finally, \(U_{1},\ldots,U_{l-1}\) are the so-called hidden layers._

Starting from any neuron \(u\in U_{\mathsf{k}}\), there is a connection to each neuron \(\hat{u}\in U_{\mathsf{k}+1}\) for \(\mathsf{k}=0,\ldots,l-1\). Such a connection illustrates that the output \(a\) of the neuron \(u\) is passed on to the neuron \(\hat{u}\). This property is the reason for the name feedforward network.

Each \(U_{0},\ldots,U_{l-1}\) contains a so-called bias neuron of the form \((0,0,1)\). It has no input, weights and a constant output value \(1\) and only transfers a constant bias in the form of the weight to each neuron of the subsequent layer.

**Remark 4.3**.: _In the following examples all neurons of the layer \(U_{\mathsf{k}}\) will have the same activation function, given by \(\sigma^{(\mathsf{k})}\) for \(\mathsf{k}=0,\ldots,l\). Here, \(D_{\mathsf{k}}:=|U_{\mathsf{k}}|-1\) for \(k=0,\ldots,l-1\) denotes the number of neurons of the \(\mathsf{k}\)-th layer (without bias neuron) and \(D_{l}:=|U_{l}|\) is the number of neurons of the output layer._

### Decomposing the domain

Before constructing the NN, we choose the domain, the decomposition and the grid on which the system (11) is solved to obtain the training values because they will influence the network size. The domain in our chosen example, given by

\[\Omega=(0,1)\times(0,1),\]

is divided into two subdomains

\[\Omega_{0}=(0,1)\times(0,0.5)\quad\text{and}\quad\Omega_{1}=(0,1)\times(0.5,1) \quad\text{(see Figure~{}\ref{fig:sub_eq_1})},\]

and the grid on which the FEM is applied is a mesh of \(32\times 32\) elements with quadratic Nedelec elements.

Hence \(32\) elements with each \(4\) degrees of freedom (dofs) are located on the interface in both subdomains. We evaluate the interface condition and the solution on each dof and use the values as the NN's input and target. Therefore the input contains \(4\cdot\dim(g_{ij})+4\cdot\dim(E_{i})=16\) values, and the output consists of \(4\cdot\dim(g_{ji})=8\) values, and we obtain \(32\) input-target pairs with one computation.

Figure 4: Visualization of the domain \(\Omega\) with the chosen decomposition.
### Neural network construction

Regarding the considerations above, we need an input layer with 16 neurons (without bias) and an output layer with 8 neurons. Furthermore, we use one hidden layer with 500 neurons (without bias). Hence, for the governing network, we have

\[U=U_{0}\dot{\cup}U_{1}\dot{\cup}U_{2}\]

with \(D_{0}=16\), \(D_{1}=500\) and \(D_{2}=8\). Our tests, presented in Section 5, revealed that this is a sufficient size for our purpose. The activation function per layer is chosen as follows:

\[\sigma^{(0)} =id\quad\text{(input layer)},\] \[\sigma^{(1)} =\frac{1}{1+e^{-x}}\quad\text{(hidden layer)},\] \[\sigma^{(2)} =id\quad\text{(output layer)},\]

where \(\sigma^{(1)}\) is known as the sigmoid function, which turned out to be the most effective since the error could be reduced more and more quickly than with other functions we tested e.g.

\[\sigma^{(1)} =\tanh(x),\] \[\sigma^{(1)} =\log\left(\frac{1}{1+e^{-x}}\right)\quad\text{(LogSigmoid)},\] \[\sigma^{(1)} =\max(0,x)+\min(0,e^{x}-1)\quad\text{(CELU)},\] \[\sigma^{(1)} =a\left(\max(0,x)+\min\left(0,b\left(e^{x}-1\right)\right) \right)\quad\text{(SELU)},\] \[\quad\text{with }a\approx 1.0507,\text{ and }b\approx 1.6733.\]

An exception represents the ReLU function, which we will discuss later in Section 5.4. Moreover, we apply separate networks \(U^{01}\) and \(U^{10}\) of the same shape for both interface conditions \(g_{01}\) and \(g_{10}\) since it turned out that they are approximated differently, fast and accurately. The resulting programming code is displayed in Figure 5.

### Training

To obtain enough training data, we vary the boundary condition \(E^{\text{inc}}\) and create training and test values to control the network during the training and avoid overfitting. The training and test sets are generated by the boundary values listed in Table 1.

Figure 5: PyTorch code of the implementation of the network construction.
Since we choose \(10\) different boundary values for the training set and \(2\) for the test set and each of them generates a set of \(32\) training/test values (one per element on the interface), we obtain all in all a set of \(32\cdot 10=320\) training values and a set of \(32\cdot 2=64\) test values for both networks. To keep the computation simple in a first set of tests, we choose a small wave number \(\omega=\frac{2\pi}{3}\), and compute the sets with the iterative DDM in \(4\) steps. Afterwards we use the results \(\left(g_{ij}^{1},E_{i}^{2}\right)\) as the input and \(g_{ji}^{3}\) as the targets to train our NNs with the application of the mean squared error as the loss function, given by

\[\text{Loss}(\mathbf{\mathsf{w}})=\frac{1}{2}\sum_{i=1}^{N}\|t^{(i)}-y(\mathbf{\mathsf{ x}}^{(i)},\mathbf{\mathsf{w}})\|^{2},\]

where \(N\) denotes the number of input-target pairs (in our case \(N=320\) for the training set and \(N=64\) for the test set), \(t^{(i)}\) is the target vector, \(y\) is the function generated by the network and hence \(y(\mathbf{\mathsf{x}}^{(i)},\mathbf{\mathsf{w}})\) denotes the output of the NN. We refer the reader to Section 5.1 for the specific realization.

As the optimizer, we use the Adam algorithm [19], which is a line search method based on the following iteration rule

\[x^{\rho+1}=x^{\rho}+\alpha^{\rho}p^{\rho},\]

where \(p^{\rho}\) is called the search direction and \(\alpha^{\rho}\) is the step size (or learning rate in case of NN) for the iteration step \(\rho\). The search direction of the Adam algorithm depends on four parameters \(\beta_{1}\), \(\beta_{2}\), \(m_{1}\) and \(m_{2}\), where \(\beta_{1}\) and \(\beta_{2}\) are fixed values in the interval \([0,1)\), and \(m_{1}\) and \(m_{2}\) are updated in each step via

\[m_{1}^{0}=m_{2}^{0}=0, m_{1}^{\rho+1}=\beta_{1}m_{1}^{\rho}+(1-\beta_{1})\cdot \nabla\text{Loss}(\mathbf{\mathsf{w}})^{\rho}\] \[\text{and} m_{2}^{\rho+1}=\beta_{2}m_{2}^{\rho}+(1-\beta_{2})\cdot\|\nabla \text{Loss}(\mathbf{\mathsf{w}})^{\rho}\|^{2}.\]

The search direction is then given by

\[p^{\rho-1}=-\widehat{m_{1}^{\rho}}\Big{/}\sqrt{\widehat{m_{2}^{\rho}}+\varepsilon}\]

with \(\widehat{m_{1}^{\rho}}=m_{1}^{\rho}/(1-(\beta_{1})^{\rho})\), \(\widehat{m_{2}^{\rho}}=m_{2}^{\rho}/(1-(\beta_{2})^{\rho})\) and \(0<\varepsilon\ll 1\).

The implementation of this training process in PyTorch is displayed in Figure 6.

\begin{table}
\begin{tabular}{|c|c|} \hline \(E^{\text{inc}}\) for the training set & \(E^{\text{inc}}\) for the test set \\ \hline \(\left(\begin{array}{c}e^{\frac{-(x-0.7)^{2}}{0.008}}\\ 0\end{array}\right)\)\(\left|\begin{array}{c}\left(\begin{array}{c}\cos(\pi^{2}y)+\sin(\pi^{2}x)i\\ \sin(\pi^{2}y)+0.5\cos(\pi^{2}x)i\end{array}\right)\end{array}\right.\) & \(\left(\begin{array}{c}e^{\frac{-(x-0.5)^{2}}{0.003}}\\ 0\end{array}\right)\) \\ \hline \(\left(\begin{array}{c}e^{\frac{-(x-0.2)^{2}}{0.002}}\\ 1\end{array}\right)\)\(\left|\begin{array}{c}\sin(\pi^{2}x)+\sin(\pi^{2}x)i\\ \sin(\pi^{2}y)+0.5\cos(\pi^{2}x)i\end{array}\right.\) & \(\left(\begin{array}{c}\cos(\pi^{2}y)+\sin(\pi^{2}x)i\\ \cos(\pi^{2}y)+0.5\cos(\pi^{2}x)i\end{array}\right)\) \\ \hline \(\left(\begin{array}{c}e^{\frac{-(x-0.7)^{2}}{0.003}}\\ 1\end{array}\right)\)\(\left|\begin{array}{c}\sin(\pi^{2}x)+\sin(\pi^{2}x)i\\ \sin(\pi^{2}x)+0.5\cos(\pi^{2}x)i\end{array}\right.\) & \\ \hline \(\left(\begin{array}{c}e^{\frac{-(x-0.8)^{2}}{0.003}}\\ \sin(\pi^{2}x)\end{array}\right)\)\(\left|\begin{array}{c}\cos(\pi^{2}y)+\sin(\pi^{2}x)i\\ \cos(\pi^{2}x)+0.5\cos(\pi^{2}x)i\end{array}\right.\) & \\ \hline \(\left(\begin{array}{c}e^{\frac{-(x-0.5)^{2}}{0.003}}\\ \cos(\pi^{2}x)\end{array}\right)\)\(\left|\begin{array}{c}\cos(\pi^{2}x)+\sin(\pi^{2}x)i\\ \cos(\pi^{2}y)+0.5\cos(\pi^{2}x)i\end{array}\right.\) & \\ \hline \end{tabular}
\end{table}
Table 1: Boundary values for generating the training set and the test set
The network \(U^{01}\) is trained with the learning rate \(10^{-5}\). The initial training error of \(3.12\) and the test error of \(5.87\) are reduced to \(1.7\cdot 10^{-4}\) and \(3\cdot 10^{-3}\) after \(29\,843\) training steps. At \(U^{10}\), the initial training error of \(0.72\) and the test error of \(1.28\) are reduced to \(3\cdot 10^{-4}\) and \(4\cdot 10^{-3}\) after \(20\,326\) steps with learning rate of \(10^{-5}\) and after further training with a learning rate of \(10^{-6}\) in \(3706\) steps, we finally achieve the training error \(2.9\cdot 10^{-4}\) and the test error \(3\cdot 10^{-3}\).

## 5 Numerical tests

In this section, we investigate several numerical experiments to demonstrate the current capacities of our approach. In addition, we highlight and analyze shortcomings and challenges.

### Comparison of new approach and classical DDM

In this first numerical example, we apply the implemented and trained NNs for the following boundary condition

\[E^{\mathrm{inc}}(x,y)=\left(\begin{array}{c}\cos\left(\pi^{2} \left(y-0.5\right)\right)+\sin\left(\pi^{2}x\right)i\\ \cos\left(\pi^{2}y\right)+0.5\sin\left(\pi^{2}x\right)i\end{array}\right),\]

and compute the first interface conditions \(g^{1}_{10}\) and \(g^{1}_{01}\) and the solutions \(E^{2}_{1}\) and \(E^{2}_{0}\) by solving (11) and (12) with the use of the parameters given in Table 2. Afterwards, these values are passed on to the networks \(U^{01}\) and \(U^{10}\). The output they return is then handled as our new interface condition, which we use to solve system (11) one more time. With that, we obtain the final solution. Moreover, we compute the same example with the DDM in 4 steps. The results that are displayed in Figure 7 show excellent agreement.

Figure 6: PyTorch code of the implementation of the network training.
### Higher wave numbers

As a second example, we increase the wave number, which leads to a more complicated problem. Therefore we repeat the same computation with \(\omega=\pi\) and leave the other parameters (especially the parameters and hyperparameters of the neural networks) unchanged. In contrast to the previous example, the results that are displayed in Figure 8 show differences. While the imaginary part is still well approximated, the real part of the NN solution differs significantly from the DDM solution and shows a discontinuity on the interface.

\begin{table}
\begin{tabular}{|c|c|c|} \hline Parameter & Definition & Value \\ \hline \hline \(\mu\) & relative magnetic permeability & 1.00 \\ \hline \(\varepsilon\) & relative electric permittivity & \(1.49^{2}\) \\ \hline \(\kappa\) & & \(\sqrt{\varepsilon}\) = 1.49 \\ \hline \(\lambda\) & wave length & 3.00 \\ \hline \(\omega\) & wave number & \(\frac{2\pi}{\lambda}+\frac{2\pi}{3.00}\) \\ \hline  & grid size & \(\frac{1}{32}\) \\ \hline \end{tabular}
\end{table}
Table 2: Parameters for the DDM

Figure 7: First example: Real part (above) and imaginary part (below) of the NN solution (left) and the DDM solution (right).
### Refined computational analysis for intermediate wave numbers

A possible reason for the mismatching results in Section 5.2 is the "problem of big wave numbers", which is very well-studied for Helmholtz-type problems [13]. The same problem also applies to Maxwell's equations [4]. To verify this conjecture and because of the very distinctive results in Section 5.1 and Section 5.2, we attempt two more computations with other wave numbers, namely \(\omega=\frac{2\pi}{2.9}\) and \(\omega=\frac{2\pi}{3.1}\). The results, that are displayed in Figures 9 and 10, in which we neglect the representation of the meshes to make the differences more visible, show that the approximation becomes inaccurate if the wave number differs slightly from the one we used for the training, regardless of whether it is larger or smaller. Therefore the bad approximation is not due to the big size of the wave number. Instead of this, it can be assumed that the NNs are specialized for the specific wave number they are trained with and "learn along" this value during the training process.

Figure 8: Second example: Real part (above) and imaginary part (below) of the NN solution (left) and the DDM solution (right).
Figure 10: Fourth example: Real part (above) and imaginary part (below) of the NN solution (left) and the DDM solution (right) with \(\omega=\frac{2\pi}{3.1}\).

Figure 9: Third example: Real part (above) and imaginary part (below) of the NN solution (left) and the DDM solution (right) with \(\omega=\frac{2\pi}{2.9}\)
### Comparison of different neural network activation functions: Sigmoid vs. ReLU

As mentioned in Section 4.3, we tested different activation functions to train the NNs before using sigmoid. One of these is the ReLU function given by

\[f(x)=\max(0,x),\]

which is implemented in the PyTorch class torch.nn.functional. This function allows a greater and faster error reduction than the others we tested, including sigmoid. In most cases, the test error of the network \(U^{01}\) can be reduced after approx. 16000 steps with a learning rate of \(10^{-5}\) and ca. 6500 steps with a learning rate of \(10^{-6}\) to \(8\cdot 10^{-4}\), which is almost a quarter compared to the final error in the training of the same NN with sigmoid as the activation function (see Section 4.4). Also, the test error of \(U^{10}\) can be reduced more quickly, namely to \(2\cdot 10^{-3}\) after ca. 3500 steps with a learning rate of \(10^{-5}\). However, we also observed that the test error grows after a short reduction phase in other cases. But in contrast, the training error continues to shrink, revealing that the training of our ReLU-networks is more susceptible to overfitting. This suspicion is strengthened when we apply the successfully trained ReLU-networks to the first example with the same procedure described in Section 5.1. The results displayed in Figure 11 show a discontinuity in the interface. This suggests that even in the lucky cases in which the test error is reduced very well, we are dealing with overfitting, and the resulting NNs cannot accurately capture the actual problem. Because of the unreliable training of the ReLU-NNs, it is reasonable to use sigmoid as the activation function instead.

Figure 11: Real part (above) and imaginary part (below) of the NN solution with the use of ReLU (left) and Sigmoid (right) as activation function.
Conclusion

In this contribution, we provided a proof of concept and feasibility study for approximating the interface operator in domain decomposition with a feedforward neural network. These concepts are applied to the time-harmonic Maxwell's equations. We carefully described the numerical framework from the algorithmic and implementation point of view. In the realization, we coupled deal.II (C++) for solving the Maxwell's equations with PyTorch for the neural network solution. Afterwards, we conducted various numerical tests that included comparing our new approach with classical domain decomposition. Then, we studied higher wave numbers in more detail. Therein, we detected difficulties, which we further investigated, revealing that the training and testing of the neural network is highly sensitive to the specific wave number. Finally, a comparison of two different neural network activation functions was undertaken. As an outlook, we plan to increase the number of subdomains to study other wave numbers further and apply the method to three-dimensional Maxwell's equations.

## Acknowledgment

This work is funded by the Deutsche Forschungsgemeinschaft (DFG) under Germany's Excellence Strategy within the Cluster of Excellence PhoenixD (EXC 2122, Project ID 390833453).
