From 077864ec1c65bdb5574704de4e870c8537fb091f Mon Sep 17 00:00:00 2001 From: Manaswinee Bezbaruah Date: Mon, 1 Nov 2021 23:57:28 -0500 Subject: [PATCH] Step-81: minor typos fixes, default parameters changed --- doc/doxygen/references.bib | 39 ++++ examples/step-81/CMakeLists.txt | 2 +- examples/step-81/doc/intro.dox | 218 +++++++++++++------ examples/step-81/doc/intro.dox-template | 252 ---------------------- examples/step-81/doc/results.dox-template | 15 -- examples/step-81/step-81.cc | 63 +++--- 6 files changed, 227 insertions(+), 362 deletions(-) delete mode 100644 examples/step-81/doc/intro.dox-template delete mode 100644 examples/step-81/doc/results.dox-template diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index 6f97c4e370..d042c7b984 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -1328,6 +1328,45 @@ url = {https://doi.org/10.1023/a:1020533003783} } +% ------------------------------------ +% Step 81 +% ------------------------------------ + +@BOOK{Schwartz1972, + AUTHOR={M. Schwartz}, + TITLE={Principles of Electrodynamics}, + Year={1972}, + PUBLISHER={McGraw-Hill Book Company}, + ADDRESS={New York}, + SERIES={International Series in Pure and Applied Physics}, +} + +@BOOK{Monk2003, + AUTHOR={P. Monk}, + TITLE={Finite Element Methods for {Maxwell's} Equations}, + YEAR={2003}, + PUBLISHER={Oxford University Press}, + SERIES={Numerical Mathematics and Scientific Computation}, +} + +@ARTICLE{Geim2004, + AUTHOR={K.S. Novoselov and A.K. Geim, S.V. Morozov and D. Jiang and Y. Zhang and S.V. Dubonos and I.V. Grigorieva and A.A. Firsov}, + TITLE={Electric Field Effect in Atomically Thin Carbon Films}, + JOURNAL={Science}, + VOLUME={306}, + PAGES={666}, + YEAR={2004} +} + +@ARTICLE{Maier2017, + author = {M. Maier and D. Margetis and M. Luskin}, + journal = {Journal of Computational Physics}, + pages = {126--145}, + title = {Dipole excitation of surface plasmon on a conducting sheet: finite element approximation and validation}, + volume = {339}, + year = {2017} +} + % ------------------------------------ % Step 82 % ------------------------------------ diff --git a/examples/step-81/CMakeLists.txt b/examples/step-81/CMakeLists.txt index c8e17bb686..d28e396a5a 100644 --- a/examples/step-81/CMakeLists.txt +++ b/examples/step-81/CMakeLists.txt @@ -23,7 +23,7 @@ SET(TARGET_SRC CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) -FIND_PACKAGE(deal.II 10.0.0 +FIND_PACKAGE(deal.II 9.2.0 HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} ) IF(NOT ${deal.II_FOUND}) diff --git a/examples/step-81/doc/intro.dox b/examples/step-81/doc/intro.dox index a4e114dc95..b9e5d19196 100644 --- a/examples/step-81/doc/intro.dox +++ b/examples/step-81/doc/intro.dox @@ -1,82 +1,155 @@ -

Introduction

- -A surface plasmon-polariton (SPP) is a slowly decaying electromagnetic wave, -excited along a conducting sheet by an electric Hertzian dipole. Suitably -rescaled time-harmonic Maxwell's equations can be used to derive a variational -form, which in turn enables a numerical observation of these SPPs by appropriate -curl-conforming finite elements. The conducting sheet is modeled as an idealized -hypersuface with an effective electric conductivity, and the weak discontinuity -for the tangential surface appears naturally in the variational form.
- -The following tutorial is a direct solver for the 2D time-harmonic Maxwell -equations describing a scattering configuration on a lower-dimensional -interface (with absorbing impedance boundary conditions). These variational -equations aim to numerically simulate a SPP on an infite sheet with with constant -isotropic conductivity embedded in two spatial dimensions. Following is a -detailed discussion on the derivation of the variational form and the appropriate -boundary conditions.
- -

Defining the Problem

- -

Time-Harmonic Maxwell's Equations

- -Consider an electromagnetic wave $(\mathbf{E},\mathbf{H})$ in a surface -$\Omega\backslash\Sigma \subset \mathbb{R}^n$, where $n=2$ or $3$. -Assume all material parameters are time-independent. - -The Maxwell's equations for this wave are, +

Introduction

+ +A surface plasmon-polariton (SPP) is a slowly decaying electromagnetic +wave, confined near a metal-air (or similar) interfaces. 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 that incorporates 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. +We discuss in particular how to set up a complex-valued (time-harmonic), +how to implement simple first-order absorbing boundary conditions and a +more sophisticated "perfectly matched layer" 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 pointers to literature. + + +

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 +@cite Schwartz1972, @cite Monk2003 : @f[ \begin{cases} --i\omega \mathbf{H} + \nabla \times \mathbf{E} = -\mathbf{M}_a,\\ -\nabla \cdot \mathbf{H} = \frac{1}{i\omega}\nabla \cdot \mathbf{M}_a,\\ -i\omega\varepsilon\mathbf{E} + \nabla\times(\mu^{-1}\mathbf{H}) = \mathbf{J}_a,\\ -\nabla\cdot(\varepsilon\mathbf{E}) = \frac{1}{i\omega}\nabla\cdot\mathbf{J}_a. +\frac{\partial}{\partial t} \mathbf{H} + \nabla \times \mathbf{E} = -\mathbf{M}_a, +\\ +\nabla \cdot \mathbf{H} = \rho, +\\ +\frac{\partial}{\partial t} (\varepsilon\mathbf{E}) - \nabla\times(\mu^{-1}\mathbf{H}) = - \mathbf{J}_a, +\\ +\nabla\cdot(\varepsilon\mathbf{E}) = \rho_m, \end{cases} @f] +where $\nabla\times\mathbf{F}(\mathbf{x})$ denotes the curl and +$\nabla\cdot\mathbf{F}(\mathbf{x})$ denotes the divergence of a vector +field $\mathbf{F}:\Omega\to\mathbb{R}^d$ and where we have set $d=2,3$. We +have introduced two (time-independent) material parameters, the +electric permittivity +$\varepsilon$ +and the +magnetic permeability +$\mu$. In addition, $\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. Both are related to their respective charge densities by a +conservation equation @cite Schwartz1972 : +@f[ +\frac{\partial}{\partial t} \rho + \nabla\cdot\mathbf{J}_a \,=\, 0, +\qquad +\frac{\partial}{\partial t} \rho_m + \nabla\cdot\mathbf{M}_a \,=\, 0. +@f] -Here, the positive parameter $\omega$ is the temporal angular frequency. -$\mathbf{J}_a$ and $\mathbf{M}_a$ are time-independent externally applied -electric-current and magnetic-current densities respectively, arising from the -time-harmonic densities $\mathcal{J}_a(x,t) = \text{Re}\{e^{-i\omega t}\mathbf{J}_a(x)\}$ -and $\mathcal{M}_a(x,t) = \text{Re}\{e^{-i\omega t}\mathbf{M}_a(x)\}$. Moreover, -$\varepsilon(x)$ and $\mu(x)$ are rank $2$ tensors representing the complex -permittivity and the relative magnetic permeability of the corresponding medium. -$\varepsilon = \varepsilon_0(x) + i\sigma(x)/\omega$, where $\varepsilon_0(x)$ is -the dielectric permittivity and $\sigma(x)$ is the surface conductivity.
-We assume some (weak) regularity of the $x$ dependent variables $(\mathbf{E}$, -$\mathbf{H})$, $(\mathbf{J}_a,\mathbf{M}_a)$, and $(\varepsilon, \mu)$ to ensure -well-posedness. - -The surface conductivity $\sigma$ gives rise to a current density, which in turn -gives rise to a jump conditions on $\Sigma$ in the tangential component (away -from the boundary) of the magnetic field. The tangential electric field is -continuous. On the idealized, oriented hypersurface $\Sigma \subset \mathbb{R}^n$, -with unit normal $\nu$ and effective surface conductivity $\sigma^{\Sigma}$, this -is modelled as, - +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 +@f[ + \mathbf{F}(\mathbf{x},t) = \text{Re}\{e^{-i\omega + t}\tilde{\mathbf{F}}(\mathbf{x})\}, +@f] +where $\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, viz., @f[ \begin{cases} -\nu \times \left[(\mu^{-1}\mathbf{H})^+ - (\mu^{-1}\mathbf{H})^-\right]|_{\Sigma} -= \sigma^{\Sigma}\left[(\nu\times \mathbf{E}\times \nu)\right]|_{\Sigma},\\ -\nu \times \left[\mathbf{E}^+ - \mathbf{E}^-\right]|_{\Sigma} = 0. +-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{cases} @f] -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}, \mathbf{H}$ approach a spherical wave uniformly in -the radial direction for points at infinity and away from the conducting sheet. - +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 @cite Geim2004. 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$ imbedded 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, viz, @f[ -\lim\limits_{|x|\to\infty} \{\mathbf{H}\times x - c^{-1}|x|\mathbf{E}\} = 0;\qquad -\lim\limits_{|x|\to\infty} \{\mathbf{E}\times x - c^{-1}|x|\mathbf{H}\} = 0;\qquad -x \not\in \Sigma + \mathbf{J}^\Sigma=\sigma^\Sigma\,\mathbf{E}_T. @f] +Here, $\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, + +@f[\text{TODO}@f] + +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 read, +@f[ +\begin{cases} +\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{cases} +@f] +Here, 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})$. -In our case, we eliminate reflection from infinity by implementing a PML and -avoid the explicit use of the last condition.

Rescaling

@@ -236,6 +309,23 @@ Find a unique $\mathbf{E} \in X(\Omega)$ such that for all $\varphi \in X(\Omega A(\mathbf{E},\varphi) = F(\varphi) @f] + +

Absorbing boundary conditions and 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}, \mathbf{H}$ approach a spherical wave uniformly in +the radial direction for points at infinity and away from the conducting sheet. + +@f[ +\lim\limits_{|x|\to\infty} \{\mathbf{H}\times x - c^{-1}|x|\mathbf{E}\} = 0;\qquad +\lim\limits_{|x|\to\infty} \{\mathbf{E}\times x - c^{-1}|x|\mathbf{H}\} = 0;\qquad +x \not\in \Sigma +@f] + +In our case, we eliminate reflection from infinity by implementing a PML and +avoid the explicit use of the last condition. +

Discretization Scheme

The variational form is discretized on a non-uniform quadrilateral mesh with diff --git a/examples/step-81/doc/intro.dox-template b/examples/step-81/doc/intro.dox-template deleted file mode 100644 index 5e67a7f8f1..0000000000 --- a/examples/step-81/doc/intro.dox-template +++ /dev/null @@ -1,252 +0,0 @@ - - -

Introduction

- -A surface plasmon-polariton (SPP) is a slowly decaying electromagnetic wave, -excited along a conducting sheet by an electric Hertzian dipole. Suitably -rescaled time-harmonic Maxwell's equations can be used to derive a variational -form, which in turn enables a numerical observation of these SPPs. - -The following tutorial is a direct solver for the 2D time-harmonic Maxwell -equations describing a scattering configuration on a lower-dimensional -interface (with absorbing impedance boundary conditions).
- -

Defining the Problem

- -

Time-Harmonic Maxwell's Equations

- -Consider an electromagnetic wave $(\mathbf{E},\mathbf{H})$ in a surface -$\Omega\backslash\Sigma \subset \mathbb{R}^n$, where $n=2 or 3. -Assume all material parameters are time-independent. -The Maxwell's equations for this wave are, - -\f{align*}{ -\begin{cases} --i\omega \mathbf{H} + \nabla \times \mathbf{E} = -\mathbf{M}_a,\\ -\nabla \cdot \mathbf{H} = \frac{1}{i\omega}\nabla \cdot \mathbf{M}_a,\\ -i\omega\varepsilon\mathbf{E} + \nabla\times(\mu^{-1}\mathbf{H}) = \mathbf{J}_a,\\ -\nabla\cdot(\varepsilon\mathbf{E}) = \frac{1}{i\omega}\nabla\cdot\mathbf{J}_a. -\f}{align*}
- -Here, the positive parameter $\omega$ is the temporal angular frequency. -$\mathbf{J}_a$ and $\mathbf{M}_a$ are time-independent externally applied -electric-current and magnetic-current densities respectively, arising from the -time-harmonic densities $\mathcal{J}_a(x,t) = \text{Re}\{e^{-i\omega t}\mathbf{J}_a(x)$ -and $\mathcal{M}_a(x,t) = \text{Re}\{e^{-i\omega t}\mathbf{M}_a(x)$. Moreover, -$\varepsilon(x)$ and $\mu(x)$ are rank 2 tensors representing the complex permittivity -and the relative magnetic permeability of the corresponding medium. -$\varepsilon = \varepsilon_0(x) + i\sigma(x)/\omega$, where $\varepsilon_0(x)$ is -the dielectric permittivity and $\sigma(x)$ is the surface conductivity.
-We assume some (weak) regularity of the $x$ dependent variables $(\mathbf{E}$, -$\mathbf{H})$, $(\mathbf{J}_a,\mathbf{M}_a)$, and $(\varepsilon, \mu)$ to ensure -well-posedness. - -The surface conductivity $\sigma$ gives rise to a current density, which in turn -gives rise to a jump conditions on $\Sigma$ in the tangential component (away -from the boundary) of the magnetic field. The tangential electric field is -continuous. On the idealized, oriented hypersurface $\Sigma \subset \mathbb{R}^n$, -with unit normal $\nu$ and effective surface conductivity $\sigma^{\Sigma}$, this -is modelled as, -\f{align*}{ -\begin{cases} -\nu \times \left[(\mu^{-1}\mathbf{H})^+ - (\mu^{-1}\mathbf{H})^-\right]|_{\Sigma} -= \sigma^{\Sigma}\left[(\nu\times \mathbf{E}\times \nu)\right]|_{\Sigma},\\ -\nu \times \left[\mathbf{E}^+ - \mathbf{E}^-\right]|_{\Sigma} = 0. -\end{cases} -\f}{align*} - -Moreover, the above equations are supplemented by the Silver-Müller radiation -condition, if the ambient (unbounded) medium is isotropic. In our case, we eliminate -reflection from infinity by implementing a PML and avoid the explicit use of this -condition. - -

Rescaling

- -We will be using a rescaled version of the Maxwell's equations where:
-1. Every length is rescaled by the free-space wavelength $2\pi k^{-1} -:= 2\pi(\omega\sqrt{\varepsilon_0\mu_0})^{-1}, where $\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$, where $J_0$ is the strength of the -prescribed dipole source at location $a$ in the $e_i$ direction in Cartesian -coordinates.
- -We introduce the rescaled variables $\mu_r$, $\varepsilon_r$, $sigma^{\Sigma}_r$, -$\hat{x}$, $\hat{\mathbf{H}}$, $\hat{\mathbf{E}}$, $\hat{\mathbf{J}}_a$, -$\hat{\mathbf{M}}_a$. - -\f{\begin{align*} -\begin{cases} -\mu_r = \frac{1}{\mu_0}\mu\\ -\varepsilon_r = \frac{1}{\varepsilon_0}\varepsilon\\ -\sigma^{Sigma}_r = \sqrt{\frac{\mu_0}{\varepsilon_0}}\sigma^{\Sigma},\\ -\hat{x} = k_0x\\ -\hat{\mathbf{H}} = \frac{k_0}{J_0}\mu^{-1}\mathbf{H}\\ -\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\\ -\hat{\mathbf{M}}_a = \frac{k_0}{\omega\mu_0 J_0}\mathbf{M}_a -\end{cases} -\end{align*}}\f - -Accordingly, our rescaled equations are: -\f{align*}{ -\begin{cases} --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}\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. -\f}{align*}
-We will omit the hat further discuss for ease of notation. - -

Variational Form

- -Let $\Omega \subset \mathbb{R}^n$ 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 $\mathbf{H}$ in -the first equation as follows: -\f{\begin{align*} -\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) -\end{align*}}\f
- -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$. Moreover, we use the subscript $T$ to denote the -tangential part of the given vector i.e. $F_T = (\nu\times F)\times\nu. We arrive -at: -\f{\begin{align*} -\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{do}x -- i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\mathbf{E})_T) \cdot -(\nabla\times(\bar{\varphi})_T)\text{d}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. -\end{align*}}\f
- -Define -\f{\begin{align*} -A(\mathbf{E},\varphi) := \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{do}x -- i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\mathbf{E})_T) \cdot -(\nabla\times(\bar{\varphi})_T)\text{d}x.\\ -F(\varphi) := 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. -\end{align*}}\f
- -Then, our rescaled weak formulation is:
-Find a unique $\mathbf{E} \in X(\Omega)$ such that for all $\varphi \in $X(\Omega)$ -\f{\begin{align*} -A(\mathbf{E},\varphi) = F(\varphi) -\end{align*}}\f
- -

Discretization Scheme

- -The variational form is discretized on a non-uniform quadrilateral mesh with -higher-order, curl-conforming Nédélec elements. This way the interface with a -weak discontinuity can be aligned with or away from the mesh, and the convergence -rate is high. Specifically, we use second-order Nédélec elements, which under our -conditions will have a convergence rate $\mathcal{O}(#dofs)$. - -Now, consider the finite element subspace $X_h(\Omega) \subset X(\Omega)$. Define -the matrices -\f{align*}{ - A_{ij} = \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 - - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_i)_T) \cdot - (\bar{\varphi}_j)_T\text{do}x - - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_i)_T) - \cdot (\nabla\times(\bar{\varphi}_j)_T)\text{d}x, -\f} -\f{align}{ - 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. -\f}
-Then under the assumption of a sufficiently refined initial mesh -the discretized variational problem is: Find a $\varphi_j \in \X_h(\Omega)$ such -that for all $\varphi_i \in X_h(\Omega)$: -\f{\begin{align*} -A_{ij} = F_i -\end{align*}}\f
- -Using a skeleton similar to step-4, we have constructed a Maxwell class and we -have used complex-valued FENedelec elements to solve our equations.
-The material coefficients such as $\mu_r^{-1}$, $J_a$, etc. are instantiated -using the Parameters class, which is also documented below. - -

Perfectly Matched Layer

-The SPP amplitude is negatively effected 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 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 essentially 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.
- -\\TODO: Insert image of the PML - -Introduce a change of coordinates -\f{\begin{align*} -x \to \bar{x} = -\begin{cases} -x + ie_r\int\limits_\rho^r s(\tau)\text{d}\tau,\;\;\;\;\;\;\; r\ge\rho -x\;\;\;\;\;\;\;\;\;\otherwise -\end{cases} -\end{align*}}\f
- -where $r = e_r \cdot x$ and $s(\tau)$ is an appropriately chosen, nonnegative -scaling function.
- -We introduce the following $2\times2$ matrices -\f{\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*}}\f
- -where -\f{\begin{align*} -d = 1 + is(r),\\ -\bar{d} = 1 + i/r \int\limits_{\rho}^{r}s(\tau)\text{d}\tau. -\end{align*}}\f
- -and $T_{e_xe_r}$ is the rotation matrix that rotates $e_r$ onto $e_x$.
- -Thus, after applying the rescaling we get the following modified parameters -\f{\begin{align*} -\bar{\mu}_r^{-1} = \frac{\mu_r^{-1}}{d},\\ -\bar{\varepsilon}_r = A^{-1} \varepsilon_r B^{-1},\\ -\bar{\sigma}^{\Sigma}_r = C^{-1} \sigma^{\Sigma}_r B^{-1}. -\end{align*}}\f - -These PML transformations are implement in our PMLParameters 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.
- -We also add an interface at $y = 0$ in our domain, in order to test the -validity of our PML.
diff --git a/examples/step-81/doc/results.dox-template b/examples/step-81/doc/results.dox-template deleted file mode 100644 index d195483004..0000000000 --- a/examples/step-81/doc/results.dox-template +++ /dev/null @@ -1,15 +0,0 @@ -

Results

- -Using the above code (without the PML) and a forcing term of a Hertzian dipole -at the center, we have generated the following solution wave. The two complex -plane solutions are followed by the two solutions in the real plane. -Furthermore, a significant resonance is observes, causing distorted images -and necessitating a PML boundary condition -TODO: add the images -The solution of the same problem with a PML of strength 2 of radii 8 and 10 -is shown below (in the same order or complex and real solutions). Clearly, -the PML significantly reduces the distortion in our solution. -TODO: add the images -Additionally, an interface is added at y = 0, and by adjusting the surface -conductivity value and the position of our dipole, we arrive at a standing wave -TODO: add the images diff --git a/examples/step-81/step-81.cc b/examples/step-81/step-81.cc index 12fb814278..9d6b7fe62c 100644 --- a/examples/step-81/step-81.cc +++ b/examples/step-81/step-81.cc @@ -66,7 +66,9 @@ #include -// @sect3{Class Template Declartions} +// @sect3{Class Template Declarations} +// We begin our actual implementation by declaring all classes with their +// data structures and methods upfront. namespace Step81 { @@ -82,24 +84,20 @@ namespace Step81 // More explanation on the use and inheritance from the ParameterAcceptor // can be found in step-60. - -/** - * - * epsilon is the Electric Permitivitty 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 inverese of the Magnetic Permiabillity 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. It is a rank 1 tensor that depends - * on the private dipole_position_, dipole_radius_, dipole_strength_, dipole_orientation_ - * variables. -*/ + // epsilon is the Electric Permitivitty 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 inverese of the Magnetic Permiabillity 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. It is a rank 1 tensor that depends + // on the private dipole_position_, dipole_radius_, dipole_strength_, dipole_orientation_ + // variables. template class Parameters : public ParameterAcceptor @@ -233,9 +231,12 @@ namespace Step81 } // @sect4{PerfectlyMatchedLayer Class} - // The PerfectlyMatchedLayer class inherits ParameterAcceptor, and it modifies our coefficients from Parameters. - // 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 PerfectlyMatchedLayer class inherits ParameterAcceptor, + // and it modifies our coefficients from Parameters. + // 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 template class PerfectlyMatchedLayer : public ParameterAcceptor @@ -284,7 +285,7 @@ namespace Step81 add_parameter("outer radius", outer_radius, "outer radius of the PML shell"); - strength = 0.; + strength = 8.; add_parameter("strength", strength, "strength of the PML"); }; @@ -362,7 +363,9 @@ namespace Step81 } - + // @sect4{Maxwell Class} + // At this point we are ready to instantiate all the major functions of + // the finite element program and also a list of variables. template class Maxwell : public ParameterAcceptor @@ -403,6 +406,9 @@ namespace Step81 // @sect4{The Constructor} + // The Constructor simply consists specifications for the mesh + // and the order of the fnite elements. These are editable through + // the .prm file. template Maxwell::Maxwell() @@ -415,7 +421,7 @@ namespace Step81 scaling = 20; add_parameter("scaling", scaling, "scale of the hypercube geometry"); - refinements = 8; + refinements = 10; add_parameter("refinements", refinements, "number of refinements of the geometry"); @@ -427,9 +433,6 @@ namespace Step81 add_parameter("quadrature order", quadrature_order, "order of the quadrature"); - - n_outputs = 2; - add_parameter("number of outputs", n_outputs, "number of output images"); } @@ -733,7 +736,7 @@ namespace Step81 { DataOut<2> data_out; data_out.attach_dof_handler(dof_handler); - data_out.add_data_vector(solution, "solution"); + data_out.add_data_vector(solution, {"real_Ex", "real_Ey", "imag_Ex", "imag_Ey"}); data_out.build_patches(); std::ofstream output("solution.vtk"); data_out.write_vtk(output); @@ -752,7 +755,7 @@ namespace Step81 } // namespace Step81 -// The following main function just calls the class step-81(), initializes the ParameterAcceptor, +// The following main function calls the class step-81(), initializes the ParameterAcceptor, // and calls the run() function. int main() -- 2.39.5