From: Manaswinee Bezbaruah Date: Tue, 22 Jun 2021 18:26:33 +0000 (-0500) Subject: Step-81: Add a skeleton X-Git-Tag: v9.4.0-rc1~136^2~24 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f4a140115644f6a8354808c5bdc9744bada414bc;p=dealii.git Step-81: Add a skeleton This commit introduces a time-harmonic Maxwell solver example step. --- diff --git a/examples/step-81/CMakeLists.txt b/examples/step-81/CMakeLists.txt new file mode 100644 index 0000000000..c8e17bb686 --- /dev/null +++ b/examples/step-81/CMakeLists.txt @@ -0,0 +1,55 @@ +## +# CMake script +## + +# Set the name of the project and target: +SET(TARGET "step-81") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) + +FIND_PACKAGE(deal.II 10.0.0 + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +# +# Are all dependencies fulfilled? +# +IF(NOT (DEAL_II_WITH_UMFPACK AND DEAL_II_WITH_COMPLEX_VALUES)) # keep in one line + MESSAGE(FATAL_ERROR " +Error! This tutorial requires a deal.II library that was configured with the following options: + DEAL_II_WITH_UMFPACK = ON + DEAL_II_WITH_COMPLEX_VALUES = ON +However, the deal.II library found at ${DEAL_II_PATH} was configured with these options: + DEAL_II_WITH_UMFPACK = ${DEAL_II_WITH_UMFPACK} + DEAL_II_WITH_COMPLEX_VALUES = ${DEAL_II_WITH_COMPLEX_VALUES} +This conflicts with the requirements." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +SET(CLEAN_UP_FILES *.log *.vtu *.pvtu) +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/examples/step-81/doc/builds-on b/examples/step-81/doc/builds-on new file mode 100644 index 0000000000..850b582a69 --- /dev/null +++ b/examples/step-81/doc/builds-on @@ -0,0 +1 @@ +step-8 diff --git a/examples/step-81/doc/intro.dox b/examples/step-81/doc/intro.dox new file mode 100644 index 0000000000..a4e114dc95 --- /dev/null +++ b/examples/step-81/doc/intro.dox @@ -0,0 +1,337 @@ + + +

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, +@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. +\end{cases} +@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, + +@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. +\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. + +@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. + +

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} +:= 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. +@f[ +\mathbf{J}_a = J_0 e_i\delta(x-a) +@f] +
+ +Accordingly, our electric permittivity and magnetic permeability are rescaled by +$\varepsilon_0$ and $\mu_0$ as follows: +@f[ +\mu_r = \frac{1}{\mu_0}\mu,\qquad +\varepsilon_r = \frac{1}{\varepsilon_0}\varepsilon. +@f] + +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: +@f[ +\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*} +@f] + +Finally, the interface conductivity is rescaled as follows: +@f[ +\sigma^{\Sigma}_r = \sqrt{\frac{\mu_0}{\varepsilon_0}}\sigma^{\Sigma}. +@f] + +Accordingly, our rescaled equations are: +@f[ +\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. +\end{cases} +@f] + +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 $\mathbf{H}$ in +the first equation as follows: +@f[ +\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) +@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$. +@f[ +\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 +- \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. +@f] + +We use the subscript $T$ to denote the tangential part of the given vector i.e. +$F_T = (\nu\times F)\times\nu$ and $[\cdot]_{\Sigma}$ to denote a jump over +$\Sigma$ i.e. $[F]_{\Sigma}(x) = \lim\limits_{s\searrow 0}(F(x+s\nu)-F(x-s\nu))$ +for $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$. +@f[ +\nu\times\mathbf{H}+\sqrt{\mu_r^{-1}\varepsilon_r}\mathbf{E}=0\qquad x\in\partial\Omega +@f] +We assume that $\mu_r^{-1}$ and $\varepsilon$ have well-defined square root. In +our numerical computation, we combine the above absorbing boundary condition +with a Perfectly Matched Layer (PML).
+ +The jump condition can be expressed as a weak discontinuity as follows: +@f[ +[\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. +@f] + +Combining, our weak form is as follows: +@f[ +\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.= +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. +@f] + +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}(curl;\Omega)$ is space of vector-valued, measurable and square +integrable functions whose (distributive) curl admits a representation by a +square integrable function. Define a Hilbert space +@f[ +X(\Omega) = \{\varphi \in \mathbf{H}(curl;\Omega)\;\;:\;\; \varphi_T|_{\Sigma} +\in L^2(\Sigma)^2,\;\varphi_T|_{\partial\Omega} \in L^2(\partial\Omega)^2\} +@f] +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 +@f[ +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{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) := 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. +@f] + +Then, our rescaled weak formulation is:
+Find a unique $\mathbf{E} \in X(\Omega)$ such that for all $\varphi \in X(\Omega)$ +@f[ +A(\mathbf{E},\varphi) = F(\varphi) +@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}(\#\text{dofs})$.
+ +Now, consider the finite element subspace $X_h(\Omega) \subset X(\Omega)$. Define +the matrices +@f[ +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{d}o_x + - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}\varphi_{i_T} + \cdot (\nabla\times \bar{\varphi}_{j_T})\;\text{d}o_x, +@f] +@f[ +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[ +A_{ij} = F_i +@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.
+ +

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.
+ + + +Introduce a change of coordinates +@f[ +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} +@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[ +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 +@f] + +where + +@f[ +d = 1 + is(r),\qquad +\bar{d} = 1 + i/r \int\limits_{\rho}^{r}s(\tau)\text{d}\tau. +@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[ +\bar{\mu}_r^{-1} = \frac{\mu_r^{-1}}{d},\qquad +\bar{\varepsilon}_r = A^{-1} \varepsilon_r B^{-1},\qquad +\bar{\sigma}^{\Sigma}_r = C^{-1} \sigma^{\Sigma}_r B^{-1}. +@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.
diff --git a/examples/step-81/doc/intro.dox-template b/examples/step-81/doc/intro.dox-template new file mode 100644 index 0000000000..5e67a7f8f1 --- /dev/null +++ b/examples/step-81/doc/intro.dox-template @@ -0,0 +1,252 @@ + + +

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/kind b/examples/step-81/doc/kind new file mode 100644 index 0000000000..e62f4e7222 --- /dev/null +++ b/examples/step-81/doc/kind @@ -0,0 +1 @@ +fluids diff --git a/examples/step-81/doc/results.dox b/examples/step-81/doc/results.dox new file mode 100644 index 0000000000..a5b174b029 --- /dev/null +++ b/examples/step-81/doc/results.dox @@ -0,0 +1,16 @@ +

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 observed, 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/doc/results.dox-template b/examples/step-81/doc/results.dox-template new file mode 100644 index 0000000000..d195483004 --- /dev/null +++ b/examples/step-81/doc/results.dox-template @@ -0,0 +1,15 @@ +

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/doc/tooltip b/examples/step-81/doc/tooltip new file mode 100644 index 0000000000..ae43afaffd --- /dev/null +++ b/examples/step-81/doc/tooltip @@ -0,0 +1 @@ +Creating a mesh. Refining it. Writing it to a file. diff --git a/examples/step-81/step-81.cc b/examples/step-81/step-81.cc new file mode 100644 index 0000000000..cf3debe59d --- /dev/null +++ b/examples/step-81/step-81.cc @@ -0,0 +1,787 @@ +/* --------------------------------------------------------------------- + * + * 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 + * Public License as published by the Free Software Foundation; either + * 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. + */ + +// @sect3{Include files} + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + + +#include +#include +#include + + + +namespace Step81 +{ + using namespace dealii; + using namespace std::complex_literals; + + + + template + class Parameters : public ParameterAcceptor + { + public: + Parameters(); + + using rank0_type = std::complex; + + using rank1_type = Tensor<1, dim, std::complex>; + + using rank2_type = Tensor<2, dim, rank0_type>; + + using curl_type = Tensor<1, dim == 2 ? 1 : dim, rank0_type>; + + public: + rank2_type epsilon(const Point &x, + types::material_id material); + + std::complex mu_inv(const Point &x, + types::material_id material); + + rank2_type sigma(const dealii::Point &x, + types::material_id left, + types::material_id right); + + rank1_type J_a(const dealii::Point &point, + types::material_id id); + + private: + rank2_type epsilon_1; + rank2_type epsilon_2; + std::complex mu_inv_1; + std::complex mu_inv_2; + rank2_type sigma_tensor; + + double dipole_radius; + Point dipole_position; + Tensor<1, dim, double> dipole_orientation; + rank0_type dipole_strength; + }; + + ////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////// + + + template + Parameters::Parameters() + : ParameterAcceptor("Parameters") + { + epsilon_1[0][0] = 1.; + epsilon_1[1][1] = 1.; + add_parameter("material 1 epsilon", + epsilon_1, + "relative permittivity of material 1"); + + epsilon_2[0][0] = 1.; + epsilon_2[1][1] = 1.; + add_parameter("material 2 epsilon", + epsilon_2, + "relative permittivity of material 2"); + + mu_inv_1 = 1.; + add_parameter("material 1 mu_inv", + mu_inv_1, + "inverse of relative permeability of material 1"); + + mu_inv_2 = 1.; + add_parameter("material 2 mu_inv", + mu_inv_2, + "inverse of relative permeability of material 2"); + + sigma_tensor[0][0] = 0.001 + 0.2i; + sigma_tensor[1][1] = 0.001 + 0.2i; + add_parameter("sigma", + sigma_tensor, + "surface conductivity between material 1 and material 2"); + + dipole_radius = 0.3; + add_parameter("dipole radius", dipole_radius, "radius of the dipole"); + + dipole_position = Point(0., 0.8); + add_parameter("dipole position", + dipole_position, + "posititon of the dipole"); + + dipole_orientation = Tensor<1, dim, double>{{0., 1.}}; + add_parameter("dipole orientation", + dipole_orientation, + "orientation of the dipole"); + + dipole_strength = 1.; + add_parameter("dipole strength", dipole_strength, "strength of the dipole"); + } + + template + typename Parameters::rank2_type + Parameters::epsilon(const Point & /*x*/, + types::material_id material) + { + return (material == 1 ? epsilon_1 : epsilon_2); + } + + template + std::complex Parameters::mu_inv(const Point & /*x*/, + types::material_id material) + { + return (material == 1 ? mu_inv_1 : mu_inv_2); + } + + template + typename Parameters::rank2_type + Parameters::sigma(const dealii::Point & /*x*/, + types::material_id left, + types::material_id right) + { + return (left == right ? rank2_type() : sigma_tensor); + } + + template + typename Parameters::rank1_type + Parameters::J_a(const dealii::Point &point, + types::material_id /*id*/) + { + 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) / + dipole_radius / dipole_radius; + J_a = dipole_strength * dipole_orientation * scale; + return J_a; + } + + + ////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////// + + + template + class PerfectlyMatchedLayer : public ParameterAcceptor + { + public: + static_assert(dim == 2, "dim == 2"); /* only works in 2D */ + + Parameters parameters; + + using rank1_type = Tensor<1, dim, std::complex>; + + using rank2_type = Tensor<2, dim, std::complex>; + + PerfectlyMatchedLayer(); + + double inner_radius; + double outer_radius; + double strength; + + std::complex d_tensor(const Point point); + + std::complex d_bar_tensor(const Point point); + + + rank2_type T_exer(std::complex d_1, + std::complex d_2, + Point point); + + rank2_type a_matrix(const Point point); + + rank2_type b_matrix(const Point point); + + rank2_type c_matrix(const Point point); + }; + + + ////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////// + + + template + PerfectlyMatchedLayer::PerfectlyMatchedLayer() + : ParameterAcceptor("PerfectlyMatchedLayer") + { + inner_radius = 12.; + add_parameter("inner radius", + inner_radius, + "inner radius of the PML shell"); + outer_radius = 15.; + add_parameter("outer radius", + outer_radius, + "outer radius of the PML shell"); + strength = 0.; + add_parameter("strength", strength, "strength of the PML"); + }; + + + template + typename std::complex + PerfectlyMatchedLayer::d_tensor(const Point point) + { + const auto radius = point.norm(); + const double s = + strength * ((radius - inner_radius) * (radius - inner_radius)) / + ((outer_radius - inner_radius) * (outer_radius - inner_radius)); + return 1 + 1.0i * s; + } + + + template + typename std::complex + PerfectlyMatchedLayer::d_bar_tensor(const Point point) + { + const auto radius = point.norm(); + const double s_bar = + strength / 3. * + ((radius - inner_radius) * (radius - inner_radius) * + (radius - inner_radius)) / + (radius * (outer_radius - inner_radius) * (outer_radius - inner_radius)); + return 1 + 1.0i * s_bar; + } + + + template + typename PerfectlyMatchedLayer::rank2_type + PerfectlyMatchedLayer::T_exer(std::complex d_1, + std::complex d_2, + Point point) + { + 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 + typename PerfectlyMatchedLayer::rank2_type + PerfectlyMatchedLayer::a_matrix(const Point point) + { + const auto d = d_tensor(point); + const auto d_bar = d_bar_tensor(point); + return invert(T_exer(d * d, d * d_bar, point)) * + T_exer(d * d, d * d_bar, point); + } + + + template + typename PerfectlyMatchedLayer::rank2_type + PerfectlyMatchedLayer::b_matrix(const Point point) + { + const auto d = d_tensor(point); + const auto d_bar = d_bar_tensor(point); + return invert(T_exer(d, d_bar, point)) * T_exer(d, d_bar, point); + } + + + template + typename PerfectlyMatchedLayer::rank2_type + PerfectlyMatchedLayer::c_matrix(const Point point) + { + const auto d = d_tensor(point); + const auto d_bar = d_bar_tensor(point); + return invert(T_exer(1. / d_bar, 1. / d, point)) * + T_exer(1. / d_bar, 1. / d, point); + } + + + ////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////// + + + template + class Maxwell : public ParameterAcceptor + { + public: + Maxwell(); + void run(); + + private: + /* run time parameters */ + double scaling; + unsigned int refinements; + unsigned int fe_order; + unsigned int quadrature_order; + bool postprocess; + unsigned int n_outputs; + + void parse_parameters_callback(); + void make_grid(); + void setup_system(); + void assemble_system(); + void solve(); + void output_results(bool postprocess); + + Parameters parameters; + PerfectlyMatchedLayer perfectly_matched_layer; + + Triangulation triangulation; + DoFHandler dof_handler; + + std::unique_ptr> fe; + + AffineConstraints constraints; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector solution; + Vector system_rhs; + }; + + + ////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////// + + + template + Maxwell::Maxwell() + : ParameterAcceptor("Maxwell") + , dof_handler(triangulation) + { + ParameterAcceptor::parse_parameters_call_back.connect( + std::bind(&Maxwell::parse_parameters_callback, this)); + + scaling = 20; + add_parameter("scaling", scaling, "scale of the hypercube geometry"); + + refinements = 8; + add_parameter("refinements", + refinements, + "number of refinements of the geometry"); + + fe_order = 0; + add_parameter("fe order", fe_order, "order of the finite element space"); + + quadrature_order = 1; + add_parameter("quadrature order", + quadrature_order, + "order of the quadrature"); + + postprocess = true; + add_parameter("postprocess", postprocess, "solution postprocessing option"); + + n_outputs = 2; + add_parameter("number of outputs", n_outputs, "number of output images"); + } + + + template + void Maxwell::parse_parameters_callback() + { + fe = std::make_unique>(FE_Nedelec(fe_order), 2); + } + + + template + void Maxwell::make_grid() + { + GridGenerator::hyper_cube(triangulation, -scaling, scaling); + triangulation.refine_global(refinements); + + 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 + void Maxwell::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); + constraints.close(); + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, + dsp, + constraints, + /* keep_constrained_dofs = */ true); + sparsity_pattern.copy_from(dsp); + system_matrix.reinit(sparsity_pattern); + } + + + template + DEAL_II_ALWAYS_INLINE inline Tensor<1, dim, std::complex> + tangential_part(const dealii::Tensor<1, dim, std::complex> &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 + void Maxwell::assemble_system() + { + QGauss quadrature_formula(quadrature_order); + QGauss face_quadrature_formula(quadrature_order); + + FEValues fe_values(*fe, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | + update_JxW_values); + FEFaceValues fe_face_values(*fe, + face_quadrature_formula, + update_values | update_gradients | + update_quadrature_points | + update_normal_vectors | + update_JxW_values); + + 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 cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + cell_matrix = 0.; + cell_rhs = 0.; + + cell->get_dof_indices(local_dof_indices); + + FEValuesViews::Vector real_part(fe_values, 0); + FEValuesViews::Vector imag_part(fe_values, dim); + const auto &quadrature_points = fe_values.get_quadrature_points(); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const Point &position = quadrature_points[q_point]; + const auto radius = position.norm(); + const auto id = cell->material_id(); + const auto inner_radius = perfectly_matched_layer.inner_radius; + + auto mu_inv = parameters.mu_inv(position, id); + auto epsilon = parameters.epsilon(position, id); + const auto J_a = parameters.J_a(position, id); + + if (radius >= inner_radius) + { + auto A = perfectly_matched_layer.a_matrix(position); + auto B = perfectly_matched_layer.b_matrix(position); + auto d = perfectly_matched_layer.d_tensor(position); + + mu_inv = mu_inv / d; + epsilon = invert(A) * epsilon * invert(B); + }; + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const auto phi_i = real_part.value(i, q_point) - + 1.0i * imag_part.value(i, q_point); + const auto curl_phi_i = real_part.curl(i, q_point) - + 1.0i * imag_part.curl(i, q_point); + + const auto rhs_value = + (1.0i * scalar_product(J_a, phi_i)) * fe_values.JxW(q_point); + cell_rhs(i) += rhs_value.real(); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + const auto phi_j = real_part.value(j, q_point) + + 1.0i * imag_part.value(j, q_point); + const auto curl_phi_j = real_part.curl(j, q_point) + + 1.0i * 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()) + { + fe_face_values.reinit(cell, face); + FEValuesViews::Vector real_part(fe_face_values, 0); + FEValuesViews::Vector imag_part(fe_face_values, dim); + + for (unsigned int q_point = 0; q_point < n_face_q_points; + ++q_point) + { + const Point position = + quadrature_points[q_point]; + const auto radius = position.norm(); + const auto id = cell->material_id(); + const auto inner_radius = + perfectly_matched_layer.inner_radius; + + + auto mu_inv = parameters.mu_inv(position, id); + auto epsilon = parameters.epsilon(position, id); + + if (radius >= inner_radius) + { + auto A = perfectly_matched_layer.a_matrix(position); + auto B = perfectly_matched_layer.b_matrix(position); + auto d = perfectly_matched_layer.d_tensor(position); + + mu_inv = mu_inv / d; + epsilon = invert(A) * epsilon * invert(B); + }; + + const auto normal = fe_face_values.normal_vector(q_point); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const auto phi_i = real_part.value(i, q_point) - + 1.0i * imag_part.value(i, q_point); + const auto phi_i_T = tangential_part(phi_i, normal); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + const auto phi_j = + real_part.value(j, q_point) + + 1.0i * 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 = + -1.0i * + scalar_product((sqrt_prod * phi_j_T), phi_i_T); + cell_matrix(i, j) += temp.real(); + } /* j */ + } /* i */ + } /* q_point */ + } + else + { + Assert(!face->at_boundary(), ExcMessage("oops!")); + 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 real_part(fe_face_values, 0); + FEValuesViews::Vector imag_part(fe_face_values, dim); + + for (unsigned int q_point = 0; q_point < n_face_q_points; + ++q_point) + { + const Point position = + quadrature_points[q_point]; + const auto radius = position.norm(); + const auto inner_radius = + perfectly_matched_layer.inner_radius; + + auto sigma = parameters.sigma(position, id1, id2); + + if (radius >= inner_radius) + { + auto B = perfectly_matched_layer.b_matrix(position); + 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 (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const auto phi_i = real_part.value(i, q_point) - + 1.0i * imag_part.value(i, q_point); + const auto phi_i_T = tangential_part(phi_i, normal); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + const auto phi_j = + real_part.value(j, q_point) + + 1.0i * imag_part.value(j, q_point); + const auto phi_j_T = tangential_part(phi_j, normal); + + const auto temp = + -1.0i * + 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 + void Maxwell::solve() + { + SparseDirectUMFPACK A_direct; + A_direct.initialize(system_matrix); + A_direct.vmult(solution, system_rhs); + } + + + template + void Maxwell::output_results(bool postprocess) + { + DataOut<2> data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "solution"); + data_out.build_patches(); + std::ofstream output("solution.vtk"); + data_out.write_vtk(output); + + + if (postprocess) + { + for (unsigned int alpha = 1; alpha <= n_outputs; alpha++) + { + std::cout << "Running step:" << alpha << std::endl; + Vector postprocessed; + postprocessed.reinit(solution); + for (unsigned int i = 0; i < dof_handler.n_dofs(); i += 2) + { + postprocessed[i] = + std::cos(-2 / n_outputs * M_PI * alpha) * solution[i] + + std::sin(-2 / n_outputs * M_PI * alpha) * solution[i + 1]; + } + data_out.add_data_vector(postprocessed, "postprocessed"); + data_out.build_patches(); + const std::string filename = + "postprocessed-" + Utilities::int_to_string(alpha) + ".vtk"; + std::ofstream output(filename); + data_out.write_vtk(output); + std::cout << "Done running step:" << alpha << std::endl; + } + } + } + + + template + void Maxwell::run() + { + make_grid(); + setup_system(); + assemble_system(); + solve(); + output_results(postprocess); + } + +} // namespace Step81 + + +int main() +{ + try + { + Step81::Maxwell<2> maxwell_2d; + dealii::ParameterAcceptor::initialize("parameters.prm"); + 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; +}