Reference documentation for deal.II version GIT d7aca55de5 2022-08-10 12:50:02+00:00
The step-81 tutorial program

This tutorial depends on step-6.

1. Introduction
2. The commented program
1. Results
2. The plain program
This program was contributed by Manaswinee Bezbaruah (Texas A&M University), and Matthias Maier (Texas A&M University).

Introduction

A surface plasmon-polariton (SPP) is a slowly decaying electromagnetic wave that is confined near a metal-air (or similar) interface. SPP structures on novel "2D" materials such as graphene, a monoatomic layer of carbon atoms arranged in a hexagonal lattice, typically have wavelengths much shorter than the wavelength of the free-space radiation. This scale separation makes SPPs on 2D materials a promising ingredient in the design of novel optical devices.

In the following, we discuss a method for observing SPPs numerically by solving a suitable electromagnetic model based on time-harmonic Maxwell's equations which incorporate jump conditions on lower-dimensional material interfaces: The conducting sheet is modeled as an idealized hypersurface with an effective electric conductivity and the weak discontinuity for the tangential surface appears naturally in the variational form.

This tutorial presents a direct solver for the time-harmonic Maxwell equations for scattering configurations with lower-dimensional interfaces. In particular, we discuss using complex values, simple first-order absorbing boundary conditions, and a more sophisticated perfectly matched layer (PML) boundary condition for electromagnetic waves.

Time-Harmonic Maxwell's Equations with interface conditions

We start the discussion with a short derivation of the governing equations and some literature references.

Derivation of time-harmonic Maxwell's equations

In two ( $$d=2$$) or three ( $$d=3$$) spatial dimensions, the time evolution of an electromagnetic wave $$(\mathbf{E},\mathbf{H})$$ that consists of an electric field component $$\mathbf{E}(t,\mathbf{x})\;:\;\mathbb{R}\times\mathbb{R}^d\to\mathbb{R}^d$$ and a magnetic field component $$\mathbf{H}(t,\mathbf{x})\;:\;\mathbb{R}\times\mathbb{R}^d\to\mathbb{R}^d$$ is described by Maxwell's equations [144], [117] :

\begin{align*} \frac{\partial}{\partial t} \mathbf{H} + \nabla \times \mathbf{E} &= -\mathbf{M}_a, \\ \nabla \cdot \mathbf{H} &= \rho_m, \\ \frac{\partial}{\partial t} (\varepsilon\mathbf{E}) - \nabla\times(\mu^{-1}\mathbf{H}) &= - \mathbf{J}_a, \\ \nabla\cdot(\varepsilon\mathbf{E}) &= \rho. \end{align*}

Here, $$\nabla\times$$ is the curl operator, $$\nabla\cdot$$ is the divergence operator, $$\varepsilon$$ is the electric permittivity, $$\mu$$ is the magnetic permeability, $$\rho$$ is the electric charge density, and $$\rho_m$$ is a corresponding (hypothetical) magnetic monopole density. $$\mathbf{J}_a$$ and $$\mathbf{M}_a$$ are the electric and magnetic flux densities which are related to their respective charge densities by the conservation equations [144]

$\frac{\partial}{\partial t} \rho + \nabla\cdot\mathbf{J}_a = 0 \text{ and } \frac{\partial}{\partial t} \rho_m + \nabla\cdot\mathbf{M}_a = 0.$

We now make the important assumption that the material parameters $$\varepsilon$$ and $$\mu$$ are time-independent and that the fields $$\mathbf{E}$$ and $$\mathbf{H}$$, the fluxes $$\mathbf{M}_a$$ and $$\mathbf{J}_a$$, as well as the densities $$\rho$$ and $$\rho_m$$ are all time-harmonic, i.e., their time evolution is completely described by

$\mathbf{F}(\mathbf{x},t) = \text{Re}\{e^{-i\omega t}\tilde{\mathbf{F}}(\mathbf{x})\},$

in which $$\omega$$ is the temporal angular frequency and $$\tilde{\mathbf{F}}(\mathbf{x})$$ is a corresponding complex-valued vector field (or density). Inserting this ansatz into Maxwell's equations, substituting the charge conservation equations and some minor algebra then yields the so-called time-harmonic Maxwell's equations:

\begin{align*} -i\omega \tilde{\mathbf{H}} + \nabla \times \tilde{\mathbf{E}} &= -\tilde{\mathbf{M}}_a, \\ \nabla \cdot \tilde{\mathbf{H}} &= \frac{1}{i\omega}\nabla \cdot \tilde{\mathbf{M}}_a, \\ i\omega\varepsilon\tilde{\mathbf{E}} + \nabla\times(\mu^{-1}\tilde{\mathbf{H}}) &= \tilde{\mathbf{J}}_a, \\ \nabla\cdot(\varepsilon\tilde{\mathbf{E}}) &= \frac{1}{i\omega}\nabla\cdot\tilde{\mathbf{J}}_a. \end{align*}

For the sake of better readability we will now drop the tilde and simply write $$\mathbf{E}(\mathbf{x})$$, $$\mathbf{H}(\mathbf{x})$$, etc., when referring to the time-harmonic fields.

Jump conditions on lower dimensional interfaces

Graphene is a two-dimensional carbon allotrope with a single atom layer that is arranged in a honeycomb lattice [126]. Due to its atomic thickness it is an example of a so-called 2D material: Compared to the other spatial dimensions (where graphene samples can reach up to several centimeters) the atomistic thickness of graphene typically ranges around 2.5 ångstrom ( $$2.5\times10^{-10}\text{m}$$). We will thus model graphene as a lower-dimensional interface $$\Sigma$$ embedded into the computational domain $$\Omega\subset\mathbb{R}^d$$. More precisely, $$\Sigma$$ is a two-dimensional sheet in three spatial dimensions, or a one-dimensional line in two spatial dimensions. The special electronic structure of graphene gives rise to a current density on the lower-dimensional interface that is modeled with an effective surface conductivity $$\sigma^\Sigma$$ obeying Ohm's Law:

$\mathbf{J}^\Sigma=\sigma^\Sigma\,\mathbf{E}_T$

in which $$\mathbf{J}^\Sigma$$ is the surface current density, $$\mathbf{E}_T$$ denotes the tangential part of the electric field $$\mathbf{E}$$, and $$\sigma^\Sigma$$ is an appropriately chosen surface conductivity that will be discussed in more detail below. The surface current density gives rise to a jump condition on $$\Sigma$$ in the tangential component of the magnetic field. This is best seen by visualizing Ampère's law:

and then taking the limit of the upper and lower part of the line integral approaching the sheet. In contrast, the tangential part of the electric field is continuous. By fixing a unit normal $$\mathbf{\nu}$$ on the hypersurface $$\Sigma$$ both jump conditions are

\begin{align*} \mathbf{\nu} \times \left[(\mu^{-1}\mathbf{H})^+ - (\mu^{-1}\mathbf{H})^-\right]|_{\Sigma} &= \sigma^{\Sigma}\left[(\mathbf{\nu}\times \mathbf{E}\times \mathbf{\nu})\right]|_{\Sigma}, \\ \mathbf{\nu} \times \left[\mathbf{E}^+ - \mathbf{E}^-\right]|_{\Sigma} &= 0. \end{align*}

The notation $$\mathbf{F}^\pm$$ indicates the limit values of the field when approaching the interface from above or below the interface: $$\mathbf{F}^\pm(\mathbf{x})=\lim_{\delta\to0,\delta>0}\mathbf{F}(\mathbf{x}\pm\delta\mathbf{\nu})$$.

Rescaling

We will be using a rescaled version of the Maxwell's equations described above. The rescaling has the following key differences:

1. Every length is rescaled by the free-space wavelength $$2\pi k^{-1} \dealcoloneq 2\pi(\omega\sqrt{\varepsilon_0\mu_0})^{-1}$$, in which $$\varepsilon_0$$ and $$\mu_0$$ denote the vacuum dielectric permittivity and magnetic permeability, respectively.
2. $$\mathbf{E}$$, $$\mathbf{H}$$, $$\mathbf{J}_a$$, $$\mathbf{M}_a$$ are all rescaled by typical electric current strength $$J_0$$, i.e., the strength of the prescribed dipole source at location $$a$$ in the $$e_i$$ direction in Cartesian coordinates (here, $$\delta$$ is the Dirac delta operator).

$\mathbf{J}_a = J_0 \mathbf{e}_i\delta(x-a)$

Accordingly, our electric permittivity and magnetic permeability are rescaled by $$\varepsilon_0$$ and $$\mu_0$$ as

$\mu_r = \frac{1}{\mu_0}\mu \text{ and } \varepsilon_r = \frac{1}{\varepsilon_0}\varepsilon.$

We use the free-space wave number $$k_0 = \omega\sqrt{\varepsilon_0\mu_0}$$ and the dipole strength, $$J_0$$ to arrive at the following rescaling of the vector fields and coordinates:

\begin{align*} \hat{x} = k_0x, &\qquad \hat{\nabla} = \frac{1}{k_0}\nabla,\\ \hat{\mathbf{H}} = \frac{k_0}{J_0}\mu^{-1}\mathbf{H},&\qquad \hat{\mathbf{E}} = \frac{k_0^2}{\omega\mu_0 J_0}\mathbf{E},\\ \hat{\mathbf{J}}_a = \frac{1}{J_0}\mathbf{J}_a,&\qquad \hat{\mathbf{M}}_a = \frac{k_0}{\omega\mu_0 J_0}\mathbf{M}_a. \end{align*}

Finally, the interface conductivity is rescaled as

$\sigma^{\Sigma}_r = \sqrt{\frac{\mu_0}{\varepsilon_0}}\sigma^{\Sigma}.$

Accordingly, our rescaled equations are

\begin{align*} -i\mu_r \hat{\mathbf{H}} + \hat{\nabla} \times \hat{\mathbf{E}} &= -\hat{\mathbf{M}}_a, \\ \hat{\nabla} \cdot (\mu_r\hat{\mathbf{H}}) &= \frac{1}{i\omega}\hat{\nabla} \cdot \hat{\mathbf{M}}_a, \\ i\varepsilon_r\hat{\mathbf{E}} + \nabla\times(\mu^{-1}\mathbf{H}) &= \mathbf{J}_a, \\ \nabla\cdot(\varepsilon\mathbf{E}) &= \frac{1}{i\omega}\hat{\nabla} \cdot\hat{\mathbf{J}}_a. \end{align*}

We will omit the hat in further discussion for ease of notation.

Variational Statement

Let $$\Omega \subset \mathbb{R}^n$$, $$(n = 2,3)$$ be a simply connected and bounded domain with Lipschitz-continuous and piecewise smooth boundary, $$\partial\Omega$$. Let $$\Sigma$$ be an oriented, Lipschitz-continuous, piecewise smooth hypersurface. Fix a normal field $$\nu$$ on $$\Sigma$$ and let $$n$$ denote the outer normal vector on $$\partial\Omega$$.

In order to arrive at the variational form, we will substitute for $$\mathbf{H}$$ in the first equation and obtain

$\nabla \times (\mu_r^{-1}\nabla\times\mathbf{E}) - \varepsilon_r \mathbf{E} = i\mathbf{J}_a - \nabla\times (\mu_r^{-1}\mathbf{M}_a).$

Now, consider a smooth test function $$\varphi$$ with complex conjugate $$\bar{\varphi}$$. Multiply both sides of the above equation by $$\bar{\varphi}$$ and integrate by parts in $$\Omega\backslash\Sigma$$.

$\int_\Omega (\mu_r^{-1}\nabla\times\mathbf{E})\cdot (\nabla\times\bar{\varphi})\;\text{d}x - \int_\Omega \varepsilon_r\mathbf{E} \cdot \bar{\varphi}\;\text{d}x - \int_\Sigma [\nu \times (\mu_r^{-1}\nabla\times\mathbf{E} + \mu^{-1}\mathbf{M}_a)]_{\Sigma}\cdot \bar{\varphi}_T\;\text{d}o_x\\ \qquad - \int_{\partial\Omega} (\nu \times (\mu_r^{-1}\nabla\times\mathbf{E} + \mu^{-1}\mathbf{M}_a)) \cdot \bar{\varphi}_T\;\text{d}o_x = i\int_\Omega \mathbf{J}_a \cdot \bar{\varphi}\;\text{d}x - \int_\Omega \mu_r^{-1}\mathbf{M}_a \cdot (\nabla \times \bar{\varphi})\;\text{d}x.$

We use the subscript $$T$$ to denote the tangential part of the given vector and $$[\cdot]_{\Sigma}$$ to denote a jump over $$\Sigma$$, i.e.,

$\mathbf{F}_T = (\mathbf{\nu}\times \mathbf{F})\times\mathbf{\nu} \text{ and } [\mathbf{F}]_{\Sigma}(\mathbf{x}) = \lim\limits_{s\searrow 0}(\mathbf{F}(\mathbf{x}+s\mathbf{\nu})-\mathbf{F}(\mathbf{x}-s\mathbf{\nu}))$

for $$\mathbf{x}\in \Sigma$$.

For the computational domain $$\Omega$$, we introduce the absorbing boundary condition at $$\partial\Omega$$, which is obtained by using a first-order approximation of the Silver-Müller radiation condition, truncated at $$\partial\Omega$$ [117].

$\nu\times\mathbf{H}+\sqrt{\mu_r^{-1}\varepsilon_r}\mathbf{E}=0\qquad x\in\partial\Omega$

We assume that $$\mu_r^{-1}$$ and $$\varepsilon$$ have well-defined square roots. In our numerical computation, we combine the above absorbing boundary condition with a PML.

The jump condition can be expressed as a weak discontinuity as follows:

$[\nu \times (\mu_r^{-1}\nabla\times\mathbf{E} + \mu^{-1}\mathbf{M}_a)]_{\Sigma} = i\sigma_r^{\Sigma}\mathbf{E}_T,\qquad \text{on }\Sigma\\ \nu \times (\mu_r^{-1}\nabla\times\mathbf{E} + \mu^{-1}\mathbf{M}_a) = -i\sqrt{\mu_r^{-1}\varepsilon_r}\mathbf{E}_T,\qquad \text{on }\partial\Omega.$

Combining, our weak form is as follows:

$\int_\Omega (\mu_r^{-1}\nabla\times\mathbf{E})\cdot (\nabla\times\bar{\varphi})\;\text{d}x - \int_\Omega \varepsilon_r\mathbf{E} \cdot \bar{\varphi}\;\text{d}x - i\int_\Sigma (\sigma_r^{\Sigma}\mathbf{E}_T) \cdot \bar{\varphi}_T\;\text{d}o_x\\ \qquad - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}\mathbf{E}_T) \cdot (\nabla\times\bar{\varphi}_T)\;\text{d}o_x.= i\int_\Omega \mathbf{J}_a \cdot \bar{\varphi}\;\text{d}x - \int_\Omega \mu_r^{-1}\mathbf{M}_a \cdot (\nabla \times \bar{\varphi})\;\text{d}x.$

Assume that $$\sigma_r^{\Sigma} \in L^{\infty}(\Sigma)^{2\times 2}$$ is matrix-valued and symmetric, and has a semidefinite real and complex part. Let $$\varepsilon_r$$ be a smooth scalar function with $$–\text{Im}(\varepsilon_r) = 0$$, or $$\text{Im}(\varepsilon_r)\ge c > 0$$ in $$\Omega$$. $$\mu_r^{-1}$$ is a smooth scalar such that $$\sqrt{\mu_r^{-1}\varepsilon_r}$$ is real valued and strictly positive in $$\partial\Omega$$.

$$\mathbf{H}(\text{curl};\Omega)$$ is space of vector-valued, measurable and square integrable functions whose weak curl admits a representation by a square integrable function. Define a Hilbert space

$X(\Omega) = \{\varphi \in \mathbf{H}(\text{curl};\Omega)\;\;:\;\; \varphi_T|_{\Sigma} \in L^2(\Sigma)^2,\;\varphi_T|_{\partial\Omega} \in L^2(\partial\Omega)^2\}$

equipped with the norm

$\|\varphi\|^2_X = \|\varphi\|^2_{L^2(\Omega)} + \|\nabla\times\varphi\|^2_{L^2(\Omega)} + \|\varphi_T\|^2_{L^2(\Sigma)} + \|\varphi_T\|^2_{L^2(\partial\Omega)}.$

Define

$A(\mathbf{E},\varphi) \dealcoloneq \int_\Omega (\mu_r^{-1}\nabla\times\mathbf{E})\cdot (\nabla\times\bar{\varphi})\;\text{d}x - \int_\Omega \varepsilon_r\mathbf{E} \cdot \bar{\varphi}\;\text{d}x - i\int_\Sigma (\sigma_r^{\Sigma}\mathbf{E}_T) \cdot \bar{\varphi}_T\;\text{d}o_x - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}\mathbf{E}_T) \cdot (\nabla\times\bar{\varphi}_T)\;\text{d}o_x.\\ F(\varphi) \dealcoloneq i\int_\Omega \mathbf{J}_a \cdot \bar{\varphi}\;\text{d}x - \int_\Omega \mu_r^{-1}\mathbf{M}_a \cdot (\nabla \times \bar{\varphi})\;\text{d}x.$

Then, our rescaled weak formulation is:

$\text{Find a unique } \mathbf{E} \in X(\Omega) \text{ such that, for all } \varphi \in X(\Omega),$

$A(\mathbf{E},\varphi) = F(\varphi).$

Absorbing boundary conditions and the perfectly matched layer

Moreover, the above equations are supplemented by the Silver-Müller radiation condition, if the ambient (unbounded) medium is isotropic. This amounts to the requirement that $$\mathbf{E}$$ and $$\mathbf{H}$$ both approach a spherical wave uniformly in the radial direction for points at infinity and away from the conducting sheet, i.e.,

$\lim\limits_{|x|\to\infty} \{\mathbf{H}\times x - c^{-1}|x|\mathbf{E}\} = 0 \text{ and } \lim\limits_{|x|\to\infty} \{\mathbf{E}\times x - c^{-1}|x|\mathbf{H}\} = 0 \text{ for } x \not\in \Sigma.$

In our case, we eliminate reflection from infinity by implementing a PML, which is described at length below, and avoid the explicit use of the last condition.

Discretization Scheme

The variational form is discretized on a non-uniform quadrilateral mesh with higher-order, curl-conforming Nédélec elements implemented by the FE_NedelecSZ class. This way the interface with a weak discontinuity is optimal, and we get optimal convergence rates.

Consider the finite element subspace $$X_h(\Omega) \subset X(\Omega)$$. Define the matrices

$A_{ij} = \int_\Omega (\mu_r^{-1}\nabla \times \varphi_j) \cdot (\nabla\times\bar{\varphi}_i)\;\text{d}x - \int_\Omega \varepsilon_r\varphi_j \cdot \bar{\varphi}_i\;\text{d}x - i\int_\Sigma (\sigma_r^{\Sigma}\varphi_{j_T}) \cdot \bar{\varphi}_{i_T}\;\text{d}o_x - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}\varphi_{j_T}) \cdot (\nabla\times \bar{\varphi}_{i_T})\;\text{d}o_x,$

$F_i = i\int_\Omega \mathbf{J}_a \cdot \bar{\varphi_i}\;\text{d}x - \int_\Omega \mu_r^{-1}\mathbf{M}_a \cdot (\nabla \times \bar{\varphi_i}) \;\text{d}x.$

Then under the assumption of a sufficiently refined initial mesh the discretized variational problem is:

$\text{Find a unique } \mathbf{E}_h = \sum_j U_j\mathbf{\varphi}_j \in X_h(\Omega) \text{ such that}$

$\sum_jA_{ij}U_j = F_i\qquad\text{for all }i.$

Perfectly Matched Layer

The SPP amplitude is negatively affected by the absorbing boundary condition and this causes the solution image to be distorted. In order to reduce the resonance and distortion in our solutions, we are implementing a Perfectly Matched Layer (PML) in the scattering configuration.

The concept of a PML was pioneered by Bérenger [41] and it is is an indispensable tool for truncating unbounded domains for wave equations and often used in the numerical approximation of scattering problems. It is essentially a thin layer with modified material parameters placed near the boundary such that all outgoing electromagnetic waves decay exponentially with no “artificial” reflection due to truncation of the domain.

Our PML is a concentric circle with modified material coefficients ( $$\varepsilon_r, \mu_r, \sigma$$). It is located in a small region near the boundary $$\partial\Omega$$ and the transformation of the material coordinates is chosen to be a function of the radial distance $$\rho$$ from the origin $$e_r$$. The normal field $$\nu$$ of $$\Sigma$$ is orthogonal to the radial direction $$e_r$$, which makes $$\mathbf{J}_a \equiv 0$$ and $$\mathbf{M}_a \equiv 0$$ within the PML.

Introduce a change of coordinates

$x \to \bar{x} = \begin{cases} x + ie_r\int\limits_\rho^r s(\tau)\text{d}\tau,\;\;\;\;\;\;\; r\ge\rho\\ x\;\;\;\;\;\;\;\;\;\text{otherwise} \end{cases}$

in which $$r = e_r \cdot x$$ and $$s(\tau)$$ is an appropriately chosen, nonnegative scaling function.

We introduce the following $$2\times2$$ matrices

\begin{align*} A &= T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}^2}, \frac{1}{d\bar{d}}\right)T_{e_xe_r} \\ B &= T_{e_xe_r}^{-1} \text{diag}\left(d,\bar{d}\right)T_{e_xe_r} \\ C &= T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}},\frac{1}{d}\right) T_{e_xe_r} \end{align*}

in which

\begin{align*} d &= 1 + is(r) \\ \bar{d} &= 1 + i/r \int\limits_{\rho}^{r}s(\tau)\text{d}\tau \end{align*}

and $$T_{e_xe_r}$$ is the rotation matrix which rotates $$e_r$$ onto $$e_x$$. Thus, after applying the rescaling, we get the following modified parameters

\begin{align*} \bar{\mu}_r^{-1} &= \frac{\mu_r^{-1}}{d}, \\ \bar{\varepsilon}_r &= A^{-1} \varepsilon_r B^{-1}, \text{ and } \\ \bar{\sigma}^{\Sigma}_r &= C^{-1} \sigma^{\Sigma}_r B^{-1}. \end{align*}

These PML transformations are implemented in our PerfectlyMatchedLayer class. After the PML is implemented, the electromagnetic wave essentially decays exponentially within the PML region near the boundary, therefore reducing reflection from the boundary of our finite domain. The decay function also depends on the strength of our PML, which can be adjusted to see more or less visible decaying in the PML region.

The commented program

Include files

The set of include files is quite standard. The most notable include is the fe/fe_nedelec_sz.h file which allows us to use the FE_NedelecSZ elements. This is an implementation of the $$H^{curl}$$ conforming Nédélec Elements that resolves the sign conflict issues that arise from parametrization (for details we refer to the documentation of the FE_NedelecSZ element).

#include <fstream>
#include <iostream>
#include <memory>

Class Template Declarations

We begin our actual implementation by declaring all classes with their data structures and methods upfront.

namespace Step81
{
using namespace dealii;
using namespace std::complex_literals;

Parameters Class

The Parameters class inherits ParameterAcceptor, and instantiates all the coefficients in our variational equations. These coefficients are passed through ParameterAcceptor and are editable through a .prm file. More explanation on the use and inheritance from the ParameterAcceptor can be found in step-60.

epsilon is the Electric Permittivity coefficient and it is a rank 2 tensor. Depending on the material, we assign the i^th diagonal element of the tensor to the material epsilon value (one of the private epsilon_1_ or epsilon_2_ variables).

mu_inv is the inverse of the Magnetic Permiability coefficient and it is a complex number.

sigma is the Surface Conductivity coefficient between material left and material right and it is a rank 2 tensor. It is only changed if we are at the interface between two materials. If we are at an interface, we assign the i^th diagonal element of the tensor to the private sigma_ value.

J_a is the strength and orientation of the dipole. As mentioned in the rescaling,

$\mathbf{J}_a = J_0 \mathbf{e}_i\delta(x-a)$

It is a rank 1 tensor that depends on the private dipole_position, dipole_radius, dipole_strength, dipole_orientation variables.

template <int dim>
class Parameters : public ParameterAcceptor
{
public:
Parameters();
using rank0_type = std::complex<double>;
using rank2_type = Tensor<2, dim, rank0_type>;
public:
rank2_type epsilon(const Point<dim> &x, types::material_id material);
std::complex<double> mu_inv(const Point<dim> & x,
types::material_id material);
rank2_type sigma(const ::Point<dim> &x,
rank1_type J_a(const ::Point<dim> &point, types::material_id id);
private:
rank2_type epsilon_1;
rank2_type epsilon_2;
std::complex<double> mu_inv_1;
std::complex<double> mu_inv_2;
rank2_type sigma_tensor;
Point<dim> dipole_position;
Tensor<1, dim> dipole_orientation;
rank0_type dipole_strength;
};
template <int dim>
Parameters<dim>::Parameters()
: ParameterAcceptor("Parameters")
{
epsilon_1[0][0] = 1.;
epsilon_1[1][1] = 1.;
epsilon_1,
"relative permittivity of material 1");
epsilon_2[0][0] = 1.;
epsilon_2[1][1] = 1.;
epsilon_2,
"relative permittivity of material 2");
mu_inv_1 = 1.;
mu_inv_1,
"inverse of relative permeability of material 1");
mu_inv_2 = 1.;
mu_inv_2,
"inverse of relative permeability of material 2");
sigma_tensor[0][0] = std::complex<double>(0.001, 0.2);
sigma_tensor[1][1] = std::complex<double>(0.001, 0.2);
sigma_tensor,
"surface conductivity between material 1 and material 2");
dipole_position = Point<dim>(0., 0.8);
add_parameter("dipole position", dipole_position, "position of the dipole");
dipole_orientation = Tensor<1, dim>{{0., 1.}};
dipole_orientation,
"orientation of the dipole");
dipole_strength = 1.;
add_parameter("dipole strength", dipole_strength, "strength of the dipole");
}
template <int dim>
typename Parameters<dim>::rank2_type
{
return (material == 1 ? epsilon_1 : epsilon_2);
}
template <int dim>
std::complex<double> Parameters<dim>::mu_inv(const Point<dim> & /*x*/,
{
return (material == 1 ? mu_inv_1 : mu_inv_2);
}
template <int dim>
typename Parameters<dim>::rank2_type
Parameters<dim>::sigma(const ::Point<dim> & /*x*/,
{
return (left == right ? rank2_type() : sigma_tensor);
}
template <int dim>
typename Parameters<dim>::rank1_type
Parameters<dim>::J_a(const ::Point<dim> &point,
{
rank1_type J_a;
const auto distance = (dipole_position - point).norm() / dipole_radius;
if (distance > 1.)
return J_a;
double scale = std::cos(distance * M_PI / 2.) *
std::cos(distance * M_PI / 2.) / (M_PI / 2. - 2. / M_PI) /
J_a = dipole_strength * dipole_orientation * scale;
return J_a;
}
Definition: point.h:111
Definition: tensor.h:503
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2087
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:190
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
unsigned int material_id
Definition: types.h:152
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)

PerfectlyMatchedLayer Class

The PerfectlyMatchedLayer class inherits ParameterAcceptor as well. It implements the transformation matrices used to modify the permittivity and permeability tensors supplied from the Parameters class. The actual transformation of the material tensors will be done in the assembly loop. The radii and the strength of the PML is specified, and the coefficients will be modified using transformation matrices within the PML region. The radii and strength of the PML are editable through a .prm file. The rotation function $$T_{exer}$$ is the same as introduced in the perfectly matched layer section of the introduction. Similarly, the matrices A, B and C are defined as follows

$A = T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}^2},\frac{1}{d\bar{d}}\right)T_{e_xe_r},\qquad B = T_{e_xe_r}^{-1} \text{diag}\left(d,\bar{d}\right)T_{e_xe_r},\qquad C = T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}},\frac{1}{d}\right) T_{e_xe_r}.\qquad$

template <int dim>
class PerfectlyMatchedLayer : public ParameterAcceptor
{
public:
static_assert(dim == 2,
"The perfectly matched layer is only implemented in 2D.");
Parameters<dim> parameters;
PerfectlyMatchedLayer();
std::complex<double> d(const Point<dim> point);
std::complex<double> d_bar(const Point<dim> point);
rank2_type rotation(std::complex<double> d_1,
std::complex<double> d_2,
rank2_type a_matrix(const Point<dim> point);
rank2_type b_matrix(const Point<dim> point);
rank2_type c_matrix(const Point<dim> point);
private:
double strength;
};
template <int dim>
PerfectlyMatchedLayer<dim>::PerfectlyMatchedLayer()
: ParameterAcceptor("PerfectlyMatchedLayer")
{
"inner radius of the PML shell");
"outer radius of the PML shell");
strength = 8.;
add_parameter("strength", strength, "strength of the PML");
}
template <int dim>
typename std::complex<double>
{
{
const double s =
return {1.0, s};
}
else
{
return 1.0;
}
}
template <int dim>
typename std::complex<double>
PerfectlyMatchedLayer<dim>::d_bar(const Point<dim> point)
{
{
const double s_bar =
strength / 3. *
return {1.0, s_bar};
}
else
{
return 1.0;
}
}
template <int dim>
typename PerfectlyMatchedLayer<dim>::rank2_type
PerfectlyMatchedLayer<dim>::rotation(std::complex<double> d_1,
std::complex<double> d_2,
{
rank2_type result;
result[0][0] = point[0] * point[0] * d_1 + point[1] * point[1] * d_2;
result[0][1] = point[0] * point[1] * (d_1 - d_2);
result[1][0] = point[0] * point[1] * (d_1 - d_2);
result[1][1] = point[1] * point[1] * d_1 + point[0] * point[0] * d_2;
return result;
}
template <int dim>
typename PerfectlyMatchedLayer<dim>::rank2_type
PerfectlyMatchedLayer<dim>::a_matrix(const Point<dim> point)
{
const auto d = this->d(point);
const auto d_bar = this->d_bar(point);
return invert(rotation(d * d, d * d_bar, point)) *
rotation(d * d, d * d_bar, point);
}
template <int dim>
typename PerfectlyMatchedLayer<dim>::rank2_type
PerfectlyMatchedLayer<dim>::b_matrix(const Point<dim> point)
{
const auto d = this->d(point);
const auto d_bar = this->d_bar(point);
return invert(rotation(d, d_bar, point)) * rotation(d, d_bar, point);
}
template <int dim>
typename PerfectlyMatchedLayer<dim>::rank2_type
PerfectlyMatchedLayer<dim>::c_matrix(const Point<dim> point)
{
const auto d = this->d(point);
const auto d_bar = this->d_bar(point);
return invert(rotation(1. / d_bar, 1. / d, point)) *
rotation(1. / d_bar, 1. / d, point);
}
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
numbers::NumberTraits< Number >::real_type norm() const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)

Maxwell Class

At this point we are ready to declare all the major building blocks of the finite element program which consists of the usual setup and assembly routines. Most of the structure has already been introduced in previous tutorial programs. The Maxwell class also holds private instances of the Parameters and PerfectlyMatchedLayers classes introduced above. The default values of these parameters are set to show us a standing wave with absorbing boundary conditions and a PML.

template <int dim>
class Maxwell : public ParameterAcceptor
{
public:
Maxwell();
void run();
private:
/* run time parameters */
double scaling;
unsigned int refinements;
unsigned int fe_order;
bool absorbing_boundary;
void parse_parameters_callback();
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results();
Parameters<dim> parameters;
PerfectlyMatchedLayer<dim> perfectly_matched_layer;
DoFHandler<dim> dof_handler;
std::unique_ptr<FiniteElement<dim>> fe;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
void run(const Iterator &begin, const typename identity< Iterator >::type &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
Definition: work_stream.h:474
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

Class Template Definitions and Implementation

The Constructor

The Constructor simply consists of default initialization a number of discretization parameters (such as the domain size, mesh refinement, and the order of finite elements and quadrature) and declaring a corresponding entry via ParameterAcceptor::add_parameter(). All of these can be modified by editing the .prm file. Absorbing boundary conditions can be controlled with the absorbing_boundary boolean. If absorbing boundary conditions are disabled we simply enforce homogeneous Dirichlet conditions on the tangential component of the electric field. In the context of time-harmonic Maxwell's equations these are also known as perfectly conducting boundary conditions.

template <int dim>
Maxwell<dim>::Maxwell()
: ParameterAcceptor("Maxwell")
, dof_handler(triangulation)
{
[&]() { parse_parameters_callback(); });
scaling = 20;
add_parameter("scaling", scaling, "scale of the hypercube geometry");
refinements = 8;
refinements,
"number of refinements of the geometry");
fe_order = 0;
add_parameter("fe order", fe_order, "order of the finite element space");
absorbing_boundary = true;
absorbing_boundary,
"use absorbing boundary conditions?");
}
template <int dim>
void Maxwell<dim>::parse_parameters_callback()
{
fe = std::make_unique<FESystem<dim>>(FE_NedelecSZ<dim>(fe_order), 2);
}
boost::signals2::signal< void()> parse_parameters_call_back

The Maxwell::make_grid() routine creates the mesh for the computational domain which in our case is a scaled square domain. Additionally, a material interface is introduced by setting the material id of the upper half ( $$y>0$$) to 1 and of the lower half ( $$y<0$$) of the computational domain to 2. We are using a block decomposition into real and imaginary matrices for the solution matrices. More details on this are available under the Results section.

template <int dim>
void Maxwell<dim>::make_grid()
{
triangulation.refine_global(refinements);
if (!absorbing_boundary)
{
for (auto &face : triangulation.active_face_iterators())
if (face->at_boundary())
face->set_boundary_id(1);
};
for (auto &cell : triangulation.active_cell_iterators())
if (cell->center()[1] > 0.)
cell->set_material_id(1);
else
cell->set_material_id(2);
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
}
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)

The Maxwell::setup_system() routine follows the usual routine of enumerating all the degrees of freedom and setting up the matrix and vector objects to hold the system data. Enumerating is done by using DoFHandler::distribute_dofs().

template <int dim>
void Maxwell<dim>::setup_system()
{
dof_handler.distribute_dofs(*fe);
std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
dof_handler,
0, /* real part */
0, /* boundary id */
constraints);
dof_handler,
dim, /* imaginary part */
0, /* boundary id */
constraints);
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
dsp,
constraints,
/* keep_constrained_dofs = */ true);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
}
void project_boundary_values_curl_conforming_l2(const DoFHandler< dim, dim > &dof_handler, const unsigned int first_vector_component, const Function< dim, number > &boundary_function, const types::boundary_id boundary_component, AffineConstraints< number > &constraints, const Mapping< dim > &mapping)
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)

This is a helper function that takes the tangential component of a tensor.

template <int dim>
tangential_part(const ::Tensor<1, dim, std::complex<double>> &tensor,
const Tensor<1, dim> & normal)
{
auto result = tensor;
result[0] = normal[1] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
result[1] = -normal[0] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
return result;
}
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:102

Assemble the stiffness matrix and the right-hand side:

\begin{align*} A_{ij} = \int_\Omega (\mu_r^{-1}\nabla \times \varphi_j) \cdot (\nabla\times\bar{\varphi}_i)\text{d}x - \int_\Omega \varepsilon_r\varphi_j \cdot \bar{\varphi}_i\text{d}x - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_j)_T) \cdot (\bar{\varphi}_i)_T\text{do}x - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_j)_T) \cdot (\nabla\times(\bar{\varphi}_i)_T)\text{d}x, \end{align*}

\begin{align} F_i = i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x - \int_\Omega \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i}) \text{d}x. \end{align}

In addition, we will be modifying the coefficients if the position of the cell is within the PML region.

template <int dim>
void Maxwell<dim>::assemble_system()
{
FEValues<dim, dim> fe_values(*fe,
FEFaceValues<dim, dim> fe_face_values(*fe,
const unsigned int dofs_per_cell = fe->dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
const unsigned int n_face_q_points = face_quadrature_formula.size();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double >> &velocity, const double factor=1.)

This is assembling the interior of the domain on the left hand side. So we are assembling

\begin{align*} \int_\Omega (\mu_r^{-1}\nabla \times \varphi_i) \cdot (\nabla\times\bar{\varphi}_j)\text{d}x - \int_\Omega \varepsilon_r\varphi_i \cdot \bar{\varphi}_j\text{d}x \end{align*}

and

\begin{align} i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x - \int_\Omega \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i}) \text{d}x. \end{align}

In doing so, we need test functions $$\varphi_i$$ and $$\varphi_j$$, and the curl of these test variables. We must be careful with the signs of the imaginary parts of these complex test variables. Moreover, we have a conditional that changes the parameters if the cell is in the PML region.

for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
FEValuesViews::Vector<dim> imag_part(fe_values, dim);
cell_rhs = 0.;
cell->get_dof_indices(local_dof_indices);
const auto id = cell->material_id();
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
auto mu_inv = parameters.mu_inv(position, id);
auto epsilon = parameters.epsilon(position, id);
const auto J_a = parameters.J_a(position, id);
const auto A = perfectly_matched_layer.a_matrix(position);
const auto B = perfectly_matched_layer.b_matrix(position);
const auto d = perfectly_matched_layer.d(position);
mu_inv = mu_inv / d;
for (const auto i : fe_values.dof_indices())
{
constexpr std::complex<double> imag{0., 1.};
const auto phi_i = real_part.value(i, q_point) -
imag * imag_part.value(i, q_point);
const auto curl_phi_i = real_part.curl(i, q_point) -
imag * imag_part.curl(i, q_point);
const auto rhs_value =
(imag * scalar_product(J_a, phi_i)) * fe_values.JxW(q_point);
cell_rhs(i) += rhs_value.real();
for (const auto j : fe_values.dof_indices())
{
const auto phi_j = real_part.value(j, q_point) +
imag * imag_part.value(j, q_point);
const auto curl_phi_j = real_part.curl(j, q_point) +
imag * imag_part.curl(j, q_point);
const auto temp =
(scalar_product(mu_inv * curl_phi_j, curl_phi_i) -
scalar_product(epsilon * phi_j, phi_i)) *
fe_values.JxW(q_point);
cell_matrix(i, j) += temp.real();
}
}
}
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
static const char A
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double >> &properties={})
Definition: generators.cc:493

Now we assemble the face and the boundary. The following loops will assemble

\begin{align*} - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_i)_T) \cdot (\bar{\varphi}_j)_T\text{do}x \end{align*}

and

\begin{align} - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_i)_T) \cdot (\nabla\times(\bar{\varphi}_j)_T)\text{d}x, \end{align}

respectively. The test variables and the PML are implemented similarly as the domain.

for (const auto &face : cell->face_iterators())
{
if (face->at_boundary())
{
const auto id = face->boundary_id();
if (id != 0)
{
fe_face_values.reinit(cell, face);
FEValuesViews::Vector<dim> imag_part(fe_face_values, dim);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
{
auto mu_inv = parameters.mu_inv(position, id);
auto epsilon = parameters.epsilon(position, id);
const auto A =
perfectly_matched_layer.a_matrix(position);
const auto B =
perfectly_matched_layer.b_matrix(position);
const auto d = perfectly_matched_layer.d(position);
mu_inv = mu_inv / d;
const auto normal =
fe_face_values.normal_vector(q_point);
for (const auto i : fe_face_values.dof_indices())
{
constexpr std::complex<double> imag{0., 1.};
const auto phi_i =
real_part.value(i, q_point) -
imag * imag_part.value(i, q_point);
const auto phi_i_T = tangential_part(phi_i, normal);
for (const auto j : fe_face_values.dof_indices())
{
const auto phi_j =
real_part.value(j, q_point) +
imag * imag_part.value(j, q_point);
const auto phi_j_T =
tangential_part(phi_j, normal) *
fe_face_values.JxW(q_point);
const auto prod = mu_inv * epsilon;
const auto sqrt_prod = prod;
const auto temp =
-imag * scalar_product((sqrt_prod * phi_j_T),
phi_i_T);
cell_matrix(i, j) += temp.real();
} /* j */
} /* i */
} /* q_point */
}
}
else
{

We are on an interior face:

const auto face_index = cell->face_iterator_to_index(face);
const auto id1 = cell->material_id();
const auto id2 = cell->neighbor(face_index)->material_id();
if (id1 == id2)
continue; /* skip this face */
fe_face_values.reinit(cell, face);
FEValuesViews::Vector<dim> imag_part(fe_face_values, dim);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
{
auto sigma = parameters.sigma(position, id1, id2);
const auto B = perfectly_matched_layer.b_matrix(position);
const auto C = perfectly_matched_layer.c_matrix(position);
sigma = invert(C) * sigma * invert(B);
const auto normal = fe_face_values.normal_vector(q_point);
for (const auto i : fe_face_values.dof_indices())
{
constexpr std::complex<double> imag{0., 1.};
const auto phi_i = real_part.value(i, q_point) -
imag * imag_part.value(i, q_point);
const auto phi_i_T = tangential_part(phi_i, normal);
for (const auto j : fe_face_values.dof_indices())
{
const auto phi_j =
real_part.value(j, q_point) +
imag * imag_part.value(j, q_point);
const auto phi_j_T = tangential_part(phi_j, normal);
const auto temp =
-imag *
scalar_product((sigma * phi_j_T), phi_i_T) *
fe_face_values.JxW(q_point);
cell_matrix(i, j) += temp.real();
} /* j */
} /* i */
} /* q_point */
}
}
constraints.distribute_local_to_global(
cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
}
}
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)

We use a direct solver from the SparseDirectUMFPACK to solve the system

template <int dim>
void Maxwell<dim>::solve()
{
A_direct.initialize(system_matrix);
A_direct.vmult(solution, system_rhs);
}
void initialize(const SparsityPattern &sparsity_pattern)
void vmult(Vector< double > &dst, const Vector< double > &src) const

The output is written into a vtk file with 4 components

template <int dim>
void Maxwell<dim>::output_results()
{
DataOut<2> data_out;
data_out.attach_dof_handler(dof_handler);
{"real_Ex", "real_Ey", "imag_Ex", "imag_Ey"});
data_out.build_patches();
std::ofstream output("solution.vtk");
data_out.write_vtk(output);
}
template <int dim>
{
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
} // namespace Step81
void write_vtk(std::ostream &out) const
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1063

The following main function calls the class step-81(), initializes the ParameterAcceptor, and calls the run() function.

int main()
{
try
{
Step81::Maxwell<2> maxwell_2d;
maxwell_2d.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
static void initialize(const std::string &filename="", const std::string &output_filename="", const ParameterHandler::OutputStyle output_style_for_output_filename=ParameterHandler::Short, ParameterHandler &prm=ParameterAcceptor::prm, const ParameterHandler::OutputStyle output_style_for_filename=ParameterHandler::DefaultStyle)

Results

The solution is written to a .vtk file with four components. These are the real and imaginary parts of the $$E_x$$ and $$E_y$$ solution waves. With the current setup, the output should read

Number of active cells: 4096
Number of degrees of freedom: 16640
Program ended with exit code: 0

Absorbing boundary conditions and the PML

The following images are the outputs for the imaginary $$E_x$$ without the interface and with the dipole centered at $$(0,0)$$. In order to remove the interface, the surface conductivity is set to 0. First, we turn off the absorbing boundary conditions and the PML. Second, we want to see the effect of the PML when absorbing boundary conditions apply. So we set absorbing boundary conditions to true and leave the PML strength to 0. Lastly, we increase the strength of the PML to 4. Change the following in the .prm file :

# use absorbing boundary conditions?
set absorbing boundary condition = false
# position of the dipole
set dipole position = 0, 0
# strength of the PML
set strength = 0
# surface conductivity between material 1 and material 2
set sigma = 0, 0; 0, 0| 0, 0; 0, 0
__global__ void set(Number *val, const Number s, const size_type N)

Following are the output images:

 Solution with no interface, Dirichlet boundary conditions and PML strength 0. Solution with no interface, absorbing boundary conditions and PML strength 0. Solution with no interface, absorbing boundary conditions and PML strength 4.

We observe that with absorbing boundary conditions and in absence of the PML, there is a lot of distortion and resonance (the real parts will not be generated without a PML). This is, as we stipulated, due to reflection from infinity. As we see, a much more coherent image is generated with an appropriate PML.

Surface Plasmon Polariton

Now, let's generate a standing wave by adding an interface at the center. In order to observe this effect, we offset the center of the dipole to $$(0, 0.8)$$ and set the surface conductivity back to $$(0.001, 0.2)$$:

# position of the dipole
set dipole position = 0, 0.8
# surface conductivity between material 1 and material 2
set sigma = 0.001, 0.2; 0, 0| 0, 0; 0.001, 0.2

Once again, we will visualize the output with absorbing boundary conditions and PML strength 0 and with absorbing boundary conditions and PML strength

1. The following tables are the imaginary part of $$E_x$$ and the real part of $$E_x$$.
 Solution with an interface, Dirichlet boundary conditions and PML strength 0. Solution with an interface, absorbing boundary conditions and PML strength 0. Solution with an interface, absorbing boundary conditions and PML strength 4.
 Solution with an interface, Dirichlet boundary conditions and PML strength 0. Solution with an interface, absorbing boundary conditions and PML strength 0. Solution with an interface, absorbing boundary conditions and PML strength 4.

The SPP is confined near the interface that we created, however without absorbing boundary conditions, we don't observe a dissipation effect. On adding the absorbing boundary conditions, we observe distortion and resonance and we still don't notice any dissipation. As expected, the PML removes the distortion and resonance. The standing wave is also dissipating and getting absorbed within the PML, and as we increase the PML strength, the standing wave will dissipate more within the PML ring.

Here are some animations to demonstrate the effect of the PML

 Solution with an interface, Dirichlet boundary conditions and PML strength 0. Solution with an interface, absorbing boundary conditions and PML strength 0. Solution with an interface, absorbing boundary conditions and PML strength 4.
 Solution with an interface, Dirichlet boundary conditions and PML strength 0. Solution with an interface, absorbing boundary conditions and PML strength 0. Solution with an interface, absorbing boundary conditions and PML strength 4.

Notes

Real and Complex Matrices

As is evident from the results, we are splitting our solution matrices into the real and the imaginary components. We started off using the $$H^{curl}$$ conforming Nédélec Elements, and we made two copies of the Finite Elements in order to represent the real and the imaginary components of our input (FE_NedelecSZ was used instead of FE_Nedelec to avoid the sign conflicts issues present in traditional Nédélec elements). In the assembly, we create two vectors of dimension $$dim$$ that assist us in extracting the real and the imaginary components of our finite elements.

Rotations and Scaling

As we see in our assembly, our finite element is rotated and scaled as follows:

const auto phi_i = real_part.value(i, q_point) - 1.0i * imag_part.value(i, q_point);

This $$\phi_i$$ variable doesn't need to be scaled in this way, we may choose any arbitrary scaling constants $$a$$ and $$b$$. If we choose this scaling, the $$\phi_j$$ must also be modified with the same scaling, as follows:

const auto phi_i = a*real_part.value(i, q_point) -
bi * imag_part.value(i, q_point);
const auto phi_j = a*real_part.value(i, q_point) +
bi * imag_part.value(i, q_point);

Moreover, the cell_rhs need not be the real part of the rhs_value. Say if we modify to take the imaginary part of the computed rhs_value, we must also modify the cell_matrix accordingly to take the imaginary part of temp. However, making these changes to both sides of the equation will not affect our solution, and we will still be able to generate the surface plasmon polariton.

cell_rhs(i) += rhs_value.imag();
cell_matrix(i) += temp.imag();

Postprocessing

We will create a video demonstrating the wave in motion, which is essentially an implementation of $$e^{-i\omega t}(Re(E) + i*Im(E))$$ as we increment time. This is done by slightly changing the output function to generate a series of .vtk files, which will represent out solution wave as we increment time. Introduce an input variable $$t$$ in the output_results() class as output_results(unsigned int t). Then change the class itself to the following:

template <int dim>
void Maxwell<dim>::output_results(unsigned int t)
{
std::cout << "Running step:" << t << std::endl;
DataOut<2> data_out;
data_out.attach_dof_handler(dof_handler);
Vector<double> postprocessed;
postprocessed.reinit(solution);
for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
{
if (i % 4 == 0)
{
postprocessed[i] = std::cos(2 * M_PI * 0.04 * t) * solution[i] -
std::sin(2 * M_PI * 0.04 * t) * solution[i + 1];
}
else if (i % 4 == 2)
{
postprocessed[i] = std::cos(2 * M_PI * 0.04 * t) * solution[i] -
std::sin(2 * M_PI * 0.04 * t) * solution[i + 1];
}
}
data_out.build_patches();
const std::string filename =
"solution-" + Utilities::int_to_string(t) + ".vtk";
std::ofstream output(filename);
data_out.write_vtk(output);
std::cout << "Done running step:" << t << std::endl;
}
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)

Finally, in the run() function, replace output_results() with

for (int t = 0; t <= 100; t++)
{
output_results(t);
}

This would generate 100 solution .vtk files, which can be opened in a group on Paraview and then can be saved as an animation. We used FFMPEG to generate gifs.

Possibilities for Extension

The example step could be extended in a number of different directions.

• The current program uses a direct solver to solve the linear system. This is efficient for two spatial dimensions where scattering problems up to a few millions degrees of freedom can be solved. In 3D, however, the increased stencil size of the Nedelec element pose a severe limiting factor on the problem size that can be computed. As an alternative, the idea to use iterative solvers can be entertained. This, however requires specialized preconditioners. For example, just using an iterative Krylov space solver (such as SolverGMRES) on above problem will requires many thousands of iterations to converge. Unfortunately, time-harmonic Maxwell's equations lack the usual notion of local smoothing properties, which renders the usual suspects, such as a geometric multigrid (see the Multigrid class), largely useless. A possible extension would be to implement an additive Schwarz preconditioner (based on domain decomposition, see for example [74]), or a sweeping preconditioner (see for example [154]).
• Another possible extension of the current program is to introduce local mesh refinement (either based on a residual estimator, or based on the dual weighted residual method, see step-14). This is in particular of interest to counter the increased computational cost caused by the scale separation between the SPP and the dipole.

The plain program

/* ---------------------------------------------------------------------
*
* Copyright (C) 2021 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Author: Manaswinee Bezbaruah, Matthias Maier, Texas A&M University, 2021.
*/
#include <fstream>
#include <iostream>
#include <memory>
namespace Step81
{
using namespace dealii;
using namespace std::complex_literals;
template <int dim>
class Parameters : public ParameterAcceptor
{
public:
Parameters();
using rank0_type = std::complex<double>;
using rank2_type = Tensor<2, dim, rank0_type>;
public:
rank2_type epsilon(const Point<dim> &x, types::material_id material);
std::complex<double> mu_inv(const Point<dim> & x,
types::material_id material);
rank2_type sigma(const ::Point<dim> &x,
rank1_type J_a(const ::Point<dim> &point, types::material_id id);
private:
rank2_type epsilon_1;
rank2_type epsilon_2;
std::complex<double> mu_inv_1;
std::complex<double> mu_inv_2;
rank2_type sigma_tensor;
Point<dim> dipole_position;
Tensor<1, dim> dipole_orientation;
rank0_type dipole_strength;
};
template <int dim>
Parameters<dim>::Parameters()
: ParameterAcceptor("Parameters")
{
epsilon_1[0][0] = 1.;
epsilon_1[1][1] = 1.;
epsilon_1,
"relative permittivity of material 1");
epsilon_2[0][0] = 1.;
epsilon_2[1][1] = 1.;
epsilon_2,
"relative permittivity of material 2");
mu_inv_1 = 1.;
mu_inv_1,
"inverse of relative permeability of material 1");
mu_inv_2 = 1.;
mu_inv_2,
"inverse of relative permeability of material 2");
sigma_tensor[0][0] = std::complex<double>(0.001, 0.2);
sigma_tensor[1][1] = std::complex<double>(0.001, 0.2);
sigma_tensor,
"surface conductivity between material 1 and material 2");
dipole_position = Point<dim>(0., 0.8);
add_parameter("dipole position", dipole_position, "position of the dipole");
dipole_orientation = Tensor<1, dim>{{0., 1.}};
dipole_orientation,
"orientation of the dipole");
dipole_strength = 1.;
add_parameter("dipole strength", dipole_strength, "strength of the dipole");
}
template <int dim>
typename Parameters<dim>::rank2_type
{
return (material == 1 ? epsilon_1 : epsilon_2);
}
template <int dim>
std::complex<double> Parameters<dim>::mu_inv(const Point<dim> & /*x*/,
{
return (material == 1 ? mu_inv_1 : mu_inv_2);
}
template <int dim>
typename Parameters<dim>::rank2_type
Parameters<dim>::sigma(const ::Point<dim> & /*x*/,
{
return (left == right ? rank2_type() : sigma_tensor);
}
template <int dim>
typename Parameters<dim>::rank1_type
Parameters<dim>::J_a(const ::Point<dim> &point,
{
rank1_type J_a;
const auto distance = (dipole_position - point).norm() / dipole_radius;
if (distance > 1.)
return J_a;
double scale = std::cos(distance * M_PI / 2.) *
std::cos(distance * M_PI / 2.) / (M_PI / 2. - 2. / M_PI) /
J_a = dipole_strength * dipole_orientation * scale;
return J_a;
}
template <int dim>
class PerfectlyMatchedLayer : public ParameterAcceptor
{
public:
static_assert(dim == 2,
"The perfectly matched layer is only implemented in 2D.");
Parameters<dim> parameters;
PerfectlyMatchedLayer();
std::complex<double> d(const Point<dim> point);
std::complex<double> d_bar(const Point<dim> point);
rank2_type rotation(std::complex<double> d_1,
std::complex<double> d_2,
rank2_type a_matrix(const Point<dim> point);
rank2_type b_matrix(const Point<dim> point);
rank2_type c_matrix(const Point<dim> point);
private:
double strength;
};
template <int dim>
PerfectlyMatchedLayer<dim>::PerfectlyMatchedLayer()
: ParameterAcceptor("PerfectlyMatchedLayer")
{
"inner radius of the PML shell");
"outer radius of the PML shell");
strength = 8.;
add_parameter("strength", strength, "strength of the PML");
}
template <int dim>
typename std::complex<double>
{
{
const double s =
return {1.0, s};
}
else
{
return 1.0;
}
}
template <int dim>
typename std::complex<double>
PerfectlyMatchedLayer<dim>::d_bar(const Point<dim> point)
{
{
const double s_bar =
strength / 3. *
return {1.0, s_bar};
}
else
{
return 1.0;
}
}
template <int dim>
typename PerfectlyMatchedLayer<dim>::rank2_type
PerfectlyMatchedLayer<dim>::rotation(std::complex<double> d_1,
std::complex<double> d_2,
{
rank2_type result;
result[0][0] = point[0] * point[0] * d_1 + point[1] * point[1] * d_2;
result[0][1] = point[0] * point[1] * (d_1 - d_2);
result[1][0] = point[0] * point[1] * (d_1 - d_2);
result[1][1] = point[1] * point[1] * d_1 + point[0] * point[0] * d_2;
return result;
}
template <int dim>
typename PerfectlyMatchedLayer<dim>::rank2_type
PerfectlyMatchedLayer<dim>::a_matrix(const Point<dim> point)
{
const auto d = this->d(point);
const auto d_bar = this->d_bar(point);
return invert(rotation(d * d, d * d_bar, point)) *
rotation(d * d, d * d_bar, point);
}
template <int dim>
typename PerfectlyMatchedLayer<dim>::rank2_type
PerfectlyMatchedLayer<dim>::b_matrix(const Point<dim> point)
{
const auto d = this->d(point);
const auto d_bar = this->d_bar(point);
return invert(rotation(d, d_bar, point)) * rotation(d, d_bar, point);
}
template <int dim>
typename PerfectlyMatchedLayer<dim>::rank2_type
PerfectlyMatchedLayer<dim>::c_matrix(const Point<dim> point)
{
const auto d = this->d(point);
const auto d_bar = this->d_bar(point);
return invert(rotation(1. / d_bar, 1. / d, point)) *
rotation(1. / d_bar, 1. / d, point);
}
template <int dim>
class Maxwell : public ParameterAcceptor
{
public:
Maxwell();
void run();
private:
/* run time parameters */
double scaling;
unsigned int refinements;
unsigned int fe_order;
bool absorbing_boundary;
void parse_parameters_callback();
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results();
Parameters<dim> parameters;
PerfectlyMatchedLayer<dim> perfectly_matched_layer;
DoFHandler<dim> dof_handler;
std::unique_ptr<FiniteElement<dim>> fe;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
template <int dim>
Maxwell<dim>::Maxwell()
: ParameterAcceptor("Maxwell")
, dof_handler(triangulation)
{
[&]() { parse_parameters_callback(); });
scaling = 20;
add_parameter("scaling", scaling, "scale of the hypercube geometry");
refinements = 8;
refinements,
"number of refinements of the geometry");
fe_order = 0;
add_parameter("fe order", fe_order, "order of the finite element space");
absorbing_boundary = true;
absorbing_boundary,
"use absorbing boundary conditions?");
}
template <int dim>
void Maxwell<dim>::parse_parameters_callback()
{
fe = std::make_unique<FESystem<dim>>(FE_NedelecSZ<dim>(fe_order), 2);
}
template <int dim>
void Maxwell<dim>::make_grid()
{
triangulation.refine_global(refinements);
if (!absorbing_boundary)
{
for (auto &face : triangulation.active_face_iterators())
if (face->at_boundary())
face->set_boundary_id(1);
};
for (auto &cell : triangulation.active_cell_iterators())
if (cell->center()[1] > 0.)
cell->set_material_id(1);
else
cell->set_material_id(2);
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
}
template <int dim>
void Maxwell<dim>::setup_system()
{
dof_handler.distribute_dofs(*fe);
std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
dof_handler,
0, /* real part */
0, /* boundary id */
constraints);
dof_handler,
dim, /* imaginary part */
0, /* boundary id */
constraints);
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
dsp,
constraints,
/* keep_constrained_dofs = */ true);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
}
template <int dim>
tangential_part(const ::Tensor<1, dim, std::complex<double>> &tensor,
const Tensor<1, dim> & normal)
{
auto result = tensor;
result[0] = normal[1] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
result[1] = -normal[0] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
return result;
}
template <int dim>
void Maxwell<dim>::assemble_system()
{
FEValues<dim, dim> fe_values(*fe,
FEFaceValues<dim, dim> fe_face_values(*fe,
const unsigned int dofs_per_cell = fe->dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
const unsigned int n_face_q_points = face_quadrature_formula.size();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
FEValuesViews::Vector<dim> imag_part(fe_values, dim);
cell_rhs = 0.;
cell->get_dof_indices(local_dof_indices);
const auto id = cell->material_id();
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
auto mu_inv = parameters.mu_inv(position, id);
auto epsilon = parameters.epsilon(position, id);
const auto J_a = parameters.J_a(position, id);
const auto A = perfectly_matched_layer.a_matrix(position);
const auto B = perfectly_matched_layer.b_matrix(position);
const auto d = perfectly_matched_layer.d(position);
mu_inv = mu_inv / d;
for (const auto i : fe_values.dof_indices())
{
constexpr std::complex<double> imag{0., 1.};
const auto phi_i = real_part.value(i, q_point) -
imag * imag_part.value(i, q_point);
const auto curl_phi_i = real_part.curl(i, q_point) -
imag * imag_part.curl(i, q_point);
const auto rhs_value =
(imag * scalar_product(J_a, phi_i)) * fe_values.JxW(q_point);
cell_rhs(i) += rhs_value.real();
for (const auto j : fe_values.dof_indices())
{
const auto phi_j = real_part.value(j, q_point) +
imag * imag_part.value(j, q_point);
const auto curl_phi_j = real_part.curl(j, q_point) +
imag * imag_part.curl(j, q_point);
const auto temp =
(scalar_product(mu_inv * curl_phi_j, curl_phi_i) -
scalar_product(epsilon * phi_j, phi_i)) *
fe_values.JxW(q_point);
cell_matrix(i, j) += temp.real();
}
}
}
for (const auto &face : cell->face_iterators())
{
if (face->at_boundary())
{
const auto id = face->boundary_id();
if (id != 0)
{
fe_face_values.reinit(cell, face);
FEValuesViews::Vector<dim> imag_part(fe_face_values, dim);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
{
auto mu_inv = parameters.mu_inv(position, id);
auto epsilon = parameters.epsilon(position, id);
const auto A =
perfectly_matched_layer.a_matrix(position);
const auto B =
perfectly_matched_layer.b_matrix(position);
const auto d = perfectly_matched_layer.d(position);
mu_inv = mu_inv / d;
const auto normal =
fe_face_values.normal_vector(q_point);
for (const auto i : fe_face_values.dof_indices())
{
constexpr std::complex<double> imag{0., 1.};
const auto phi_i =
real_part.value(i, q_point) -
imag * imag_part.value(i, q_point);
const auto phi_i_T = tangential_part(phi_i, normal);
for (const auto j : fe_face_values.dof_indices())
{
const auto phi_j =
real_part.value(j, q_point) +
imag * imag_part.value(j, q_point);
const auto phi_j_T =
tangential_part(phi_j, normal) *
fe_face_values.JxW(q_point);
const auto prod = mu_inv * epsilon;
const auto sqrt_prod = prod;
const auto temp =
-imag * scalar_product((sqrt_prod * phi_j_T),
phi_i_T);
cell_matrix(i, j) += temp.real();
} /* j */
} /* i */
} /* q_point */
}
}
else
{
const auto face_index = cell->face_iterator_to_index(face);
const auto id1 = cell->material_id();
const auto id2 = cell->neighbor(face_index)->material_id();
if (id1 == id2)
continue; /* skip this face */
fe_face_values.reinit(cell, face);
FEValuesViews::Vector<dim> imag_part(fe_face_values, dim);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
{
auto sigma = parameters.sigma(position, id1, id2);
const auto B = perfectly_matched_layer.b_matrix(position);
const auto C = perfectly_matched_layer.c_matrix(position);
sigma = invert(C) * sigma * invert(B);
const auto normal = fe_face_values.normal_vector(q_point);
for (const auto i : fe_face_values.dof_indices())
{
constexpr std::complex<double> imag{0., 1.};
const auto phi_i = real_part.value(i, q_point) -
imag * imag_part.value(i, q_point);
const auto phi_i_T = tangential_part(phi_i, normal);
for (const auto j : fe_face_values.dof_indices())
{
const auto phi_j =
real_part.value(j, q_point) +
imag * imag_part.value(j, q_point);
const auto phi_j_T = tangential_part(phi_j, normal);
const auto temp =
-imag *
scalar_product((sigma * phi_j_T), phi_i_T) *
fe_face_values.JxW(q_point);
cell_matrix(i, j) += temp.real();
} /* j */
} /* i */
} /* q_point */
}
}
constraints.distribute_local_to_global(
cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
}
}
template <int dim>
void Maxwell<dim>::solve()
{
A_direct.initialize(system_matrix);
A_direct.vmult(solution, system_rhs);
}
template <int dim>
void Maxwell<dim>::output_results()
{
DataOut<2> data_out;
data_out.attach_dof_handler(dof_handler);
{"real_Ex", "real_Ey", "imag_Ex", "imag_Ey"});
data_out.build_patches();
std::ofstream output("solution.vtk");
data_out.write_vtk(output);
}
template <int dim>
{
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
} // namespace Step81
int main()
{
try
{
Step81::Maxwell<2> maxwell_2d;
maxwell_2d.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}