Medusa: Coordinate Free Mehless Method implementation - User contributions [en]

Computational electromagnetics

2020-09-27T15:21:08Z

BlazS: /* Streched Coordinate Perfectly Matched Layer */

== Case studies ==
The following pages describe the basics of computational electromagnetics, starting with relevant derivations and the basics of classical electromagnetism. The subpages include case studies with analytical solutions if they exist and numerical solutions.
* [[Triple dielectric step in 1D]]
* [[Scattering from an infinite cylinder]]
* [[Point source near an anisotropic lens]]

== Classical electromagnetism ==
=== Maxwell's equations in matter ===
Classical electrodynamics is historically one of the most eminent fields of physics as an extension of classical mechanics, it is very successful in explaining a plethora of phenomena. The dynamics of electric and magnetic fields are described with Maxwell's equations. As we will be studying the interaction of electromagnetic waves with different objects, we need Maxwell's equations in matter
<math>
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\begin{align}
&\nabla \times \b{E}(\b{r}, t) = - \dpar{\b{B}(\b{r}, t)}{t}, \label{eq:TFaraday} \\
&\nabla \times \b{H}(\b{r}, t) = \b{j}(\b{r}, t) + \dpar{\b{D}(\b{r}, t)}{t}, \label{eq:TMaxwell-Ampere} \\
&\nabla \cdot \b{D}(\b{r}, t) = \rho(\b{r}, t), \label{eq:TGaussE} \\
&\nabla \cdot \b{B}(\b{r}, t) = 0. \label{eq:TGaussM}
\end{align}
</math>

The system of equations contains four fields. The electric field strength $\b E$ and density $\b D$ and the magnetic field strength $\b H$ and density $\b B$. The four fields are accompanied by the current density $\b j$ and the charge density $\rho$. For a full description of electromagnetic phenomena, we need to provide another two constitutive equations, that relate the strength and density of the fields
<math>
\newcommand{\eps}{\varepsilon}
\begin{align}
\b{B} &= \mu_0 \mu \b{H}, \label{eq:constM} \\
\b{D} &= \eps_0 \eps \b{E} \label{eq:constE}
\end{align}
</math>
where $\varepsilon_0$ and $\mu_0$ are vacuum permittivity and permeability respectively. The dielectric function $\varepsilon$ and magnetic permeability $\mu$ are in general dependant on both $\b E$ and $\b B$ as well as the frequency $\omega$ and can be second order tensors in anisotropic materials. Equations \eqref{eq:constM} and \eqref{eq:constE} already assume linear material properties, generally the polarisation $\b P$ and magnetisation are defined as power series expansions of the electric field density and magnetic field density, as
<math>
\begin{align}
\label{eq:PMexpans}
\b{P} &= \chi_E \b{D} + \mathcal{O}(D^2), \\
\b{M} &= \chi_M \b{H} + \mathcal{O}(H^2).
\end{align}
</math>
The material linearity assumption holds well for small external fields, meaning small $\b D$ and $\b H$. The treatment of nonlinear terms falls within the field of nonlinear optics and is not relevant for our discussion here.

=== Electromagnetic waves ===
<math>
\def\doubleunderline#1{\underline{\underline{#1}}}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\newcommand{\eps}{\varepsilon}
</math>
When studying light and related phenomena, we are usually interested in the behavior of electromagnetic waves, where the electric and magnetic fields oscillate with constant frequency and propagate in a single direction. The wave formulation helps with understanding optical phenomena such as diffraction and reflection and offers easily graspable quantities such as wave amplitude and wavelength. Wave equations for $\b E$ and $\b H$ can be derived from Maxwell's equations for linear and homogeneous materials. Let's start by treating a more general case, where $\doubleunderline \varepsilon$ is a constant second-order tensor and $\mu$ is a scalar quantity. From Maxwell's equations and constitutive relations we obtain
\begin{equation}
\nabla \times \left[ \doubleunderline{\varepsilon}^{-1} \nabla \times \b{H} \right] = \nabla \times \left[ \doubleunderline{\varepsilon}^{-1} \b j \right] - \mu \mu_0 \varepsilon_0 \ddpar{\b H}{t}.
\end{equation}
and
\begin{equation}
\nabla \times \nabla \times \b{E} = \nabla (\nabla \cdot \b E) - \nabla^2 \b E = -\mu\mu_0\dpar{\b j}{t} - \eps_0 \doubleunderline{\varepsilon} \mu \mu_0 \ddpar{\b E}{t}
\end{equation}
which simplify to wave equations for homogenous (linear) materials and no external currents and charges
\begin{equation}
\nabla^2 \b E = \eps\eps_0\mu\mu_0\ddpar{\b E}{t}, \qquad \text{ or } \qquad \nabla^2 \b H = \eps\eps_0\mu\mu_0\ddpar{\b H}{t}.
\end{equation}
One possible way of solving electromagnetic problems is simply to solve the above wave equations numerically. Another approach is to solve Maxwell's equations in the frequency domain. We start by Fourier expansion of the fields $\b E$, $\b D$, $\b B$ and $\b H$ analogous to
\begin{equation}
\label{eq:Fdecomp}
\b{E}(\b{r}, t) = \int\frac{d\omega}{2\pi} \b{E}(\b{r}, \omega) e^{-i\omega t},
\end{equation}
in the frequency domainall the fields become complex. The ratio between the real and complex component of the fields is the phase difference of the material response to external fields. The Fourier expansion leads to the harmonic form of Maxwell's equations
\begin{equation}
\label{eq:FFaraday}
\nabla \times \b{E}(\b{r}, \omega) = - i\omega \b{B} (\b{r}, \omega),
\end{equation}
\begin{equation}
\label{eq:FMaxwell-Ampere}
\nabla \times \b{H}(\b{r}, \omega) = i\omega \b{D}(\b{r}, \omega).
\end{equation}
The constitutive relations are unchanged. If we assume that there are currents we obtain a wave equation for either $\b D$ or $\b H$,
\begin{equation}
\label{eq:frekaniwave}
\nabla \times \left[ \doubleunderline{\eps}^{-1} \nabla \times \b{H} \right] = \omega^2 \mu_0 \eps_0 \mu \b{H}.
\end{equation}
The above form allows for a calculation of a steady state response to an incident harmonic plane wave, meaning that we can study scattering of any incident wave that can be represented by a plane wave expansion. The time dependent solution is obtained by an inverse Fourier transform. In empty space \eqref{eq:frekaniwave} simplifies to the Helmholtz equation
\begin{equation}
\label{eq:frekemptywave}
\nabla^2 \b{H} = - \omega^2 \mu_0 \eps_0 \b{H} = - k^2 \b{H}.
\end{equation}
Wavenumber $k$ in a vacuum is defined as $k = \omega \sqrt{\mu_0 \epsilon_0} = \frac{2 \pi}{\lambda}$. One of the major advantages of the frequency domain approach over the time domain approach is that the frequency dependence of the dielectric function does not need to be considered, as one simply uses the value of the dielectric function at a certain frequency.

===Boundary conditions for Maxwell's equations===
For modeling purposes the treatment of electromagnetic fields on boundaries between materials with different properties, such as $\varepsilon$ and $\mu$, is crucial. The boundary conditions for Maxwell's equations between two linear and homogenous mediums are as follows
\begin{equation}
\begin{array}{rl}
D_{1 n}-D_{2 n}=\sigma, \qquad & E_{1 t}-E_{2 t}=0 \\
B_{1 n}-B_{2 n}=0, \qquad & H_{1 t}-H_{2 t}=K.\end{array}
\end{equation}
Where we have fields $\b{E}_1, \b{D}_1, \b{B}_1, \b{H}_1$ in medium $1$ and fields $\b{E}_2, \b{D}_2, \b{B}_2, \b{H}_2$ in medium $2$. Subscript $t$ and $n$ denote the tangential and normal component respectively, $\sigma$ is the surface charge density and $\b K$ is the surface current density.

===Frequency dependence of the dielectric function===
As mentioned above the dielectric function $\varepsilon$ is a function of frequency. Some analytical properties of the dielectric function $\varepsilon(\omega)$ are well known, due to their relevance in optics, one example being the Kramers-Kronig relations. The dielectric function is in general, based on a nonlocal response function $\chi$, defined as
\begin{equation}
\varepsilon(\omega) = 1 + \int_{0}^{\infty} \chi(\tau) e^{i \omega \tau} \text{d}\tau.
\end{equation}
For a more detailed analysis of the analytical properties of the function, refer to any standard optics or advanced electromagnetics textbook. When it comes to numerical simulations, models of the dielectric function are of more interest than its analytical properties.
====Dielectric materials====
The simplest model for dielectrics is based on the equation of motion for the bound charges
\begin{equation}
m \frac{\mathrm{d}^{2} \mathbf{r}}{\mathrm{d} t^{2}}=-m \gamma \dot{\mathbf{r}}-m \omega_{0}^{2} \mathbf{r}+e \mathbf{E}(t)
\end{equation}
the terms on the right side describe the dissipation of energy, oscillations and the driving force of the external electric field. A short calculation leads to an expression of the frequency dependence of the dielectric function
\begin{equation}
\varepsilon(\omega)=1+ \frac{\omega_p}{\left(\omega_{0}^{2}-\omega^{2}\right)-i \gamma \omega}.
\end{equation}
The above expression can be simplified for high and low frequencies leading to the so-called Debye and Plasmonic relaxation models. The Lorentz model applies for frequencies where no terms in the denominator can be neglected. Typical Debye and Lorentz relaxations are shown in the plots below.

[[File:Debye.png|400px]][[File:Lorentz.png|400px]]

In general, the dielectric function of a material is modeled as a sum of many Debye and Lorentz relaxations, the parameters of which are determined by microscopic properties of the material itself, an extremely interesting case study is the dielectric function of water.
====Conductors====
When dealing with conductors we use the Drude model of conduction
\begin{equation}
m \frac{d \mathbf{v}(t)}{d t}=-m \gamma \mathbf{v}(t)+e \mathbf{E}(t)
\end{equation}
from which we quickly derive the equation for conductivity
\begin{equation}
\sigma(\omega)=\frac{n e^{2}}{m} \frac{1}{\gamma-i \omega},
\end{equation}
which we can relate to the dielectric function
\begin{equation}
\varepsilon(\omega)=1-\frac{\sigma(\omega)}{i \varepsilon_{0} \omega}.
\end{equation}
==Modeling in the time domain==
The single most frequent type of calculation used in computational electromagnetics is the simulation of electromagnetic fields in the time domain. The calculation propagates $\b E(\b r, t)$ and $\b H(\b r, t)$ in time, usually starting with all fields set to zero and some time dependent current source $\b J(\b r, t)$. The main advantage of time domain approaches is the ability to study transient phenomena while its main drawback is the study resonant problems.
===Finite differences in time domain (FDTD)===
==== Formulation ====
Finite differences in time domain were suggested in 1966 by Yee. The method uses finite differences to discretize Maxwell's equations and propagates the fields in time by alternating between updating the electric field strength and magnetic field strength. One of its defining attributes is the use of a staggered grid to represent $\b H$ and $\b E$ nodes. The staggered grid simplifies the calculation of curl expressions and helps with numerical stability and treating the boundaries of materials.

[[File:FDTD2.png|400px|center]]

This type of solution procedure provides an alternative to solving the wave equation for either $\b H$ or $\b E$.

==== Modeling sources ====
Since we are starting with all fields set to zero, we need to somehow introduce power into the model. This is usually done with current sources. The current signal used depends on the type of calculation. If we are interested in a steady-state solution we would use a sinusoidal source, that is slowly turned on and simulated the device until a steady state is reached. On the other hand, if we are interested in calculating the response of a device to many frequencies, we would use a Gaussian pulse and simulate the device until all the energy has left the system.

Taflove in his book [1] differentiates between five types of sources. The easiest to implement are so-called ''hard sources'', where one simply sets the value of a field in a node to the source value e.g.
\begin{equation}
\b E_z(i, t) = \sin 2\pi ft.
\end{equation}
The problem is that the node now acts as a perfect dielectric and incoming waves reflect from the point and do not pass through as they should. On the other hand one can use a ''soft source'' where energy is inserted into the domain via electric and magnetic (non-physical) current
\begin{equation}
\b E_z(i, t) = E_z(i, t) + J(t) \qquad \text{in} \qquad \b H_x(i, t) = H_x(i, t) + M(t).
\end{equation}
Waves no longer reflect from the source, however, the source still radiates in all directions, which is sometimes not desired. This issue is solved with TFSF (total field, scattered field) sources, which are based on the linearity of Maxwell's equations and allow to construct sources that inject the energy only in the desired direction. This is accomplished by decomposing the total fields into the scattered and incident
\begin{equation}
\b{E} = \b{E}^i + \b{E}^s \text{, } \qquad \b{H} = \b{H}^i + \b{H}^s.
\end{equation}
and then dividing the computational domain $\Omega$ into two subdomains. In one of the subdomains, we calculate only using the scattered part of the fields and in the other, we calculate with the total field. With careful treatment of the boundary between the subdomains one can construct an arbitrary incident wave. Alongside the already mentioned source types, ''waveguide'' sources and ''total scattered field'' sources. For an in-depth discussion refer to Taflove [1].

==== Transmittance and Reflectivity ====
Transmittance and reflectivity are dependent on the frequency of the incoming wave. One of the main strengths of time domain calculations is that one can calculate both coefficients for a wide frequency range with a single simulation. This is achieved by using a source that contains a wide range of frequencies i.e. a short Gaussian pulse.

Let us illustrate how to calculate $R$ and $T$ in one dimension. The calculation is similar in higher dimensions. We make two assumptions, firstly the source is directed directly into the device. This means that we will have only the reflected field on one side of the device and only the transmitted field on the other side of the device. Secondly, we assume that nothing reflects from the edges of the domain. We record the signals on both edges of the domain and label them $E_r(t)$ for reflected and $E_t(t)$ for transmitted. It is important to simulate for a long enough time that there is no energy left in the device. Otherwise, the calculation of the coefficients has not converged yet. The Fourier transforms of the signals at the edges, along with the transform of the source, define $R$ and $T$
\begin{equation}
R(\omega) = \left| \frac{E_o(\omega)}{S(\omega)} \right| \qquad \text{and} \qquad T(\omega)= \left| \frac{E_p(\omega)}{S(\omega)} \right|.
\end{equation}
Since we only know the field values at discrete time intervals the transformation is discrete
\begin{equation}
F(\omega) \cong \Delta t \sum_{m=1}^{M}\left(e^{-j \omega \Delta t}\right)^{m} \cdot f(m \Delta t),
\end{equation}
in practice, so that it is not needed to store the signals on both edges which may take up considerable computer memory, the sum can be calculated after each step in the simulation to save space.

<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}
</math>

==Modeling in the frequency domain==
In the frequency domain, two types of calculations are most common. Within the field of photonic crystals, bandstructures of periodic materials are sought by solving a generalized eigenvalue problem $Ax = \omega^2Bx$. In engineering applications, a response of a structure to a periodical source $\b J(x) e^{-i\omega t}$ is more often calculated. We will focus on the second type of calculation.

==== Formulation ====
We need to solve the equation that can be derived by combining Maxwell's equations and constitutive relations
\begin{equation}
\left[(\nabla \times \nabla \times)-\frac{\omega^{2}}{c^{2}} \varepsilon(\mathbf{r})\right] \mathbf{E}(\mathbf{r})=i \omega \mu_{0} \mathbf{J}(\mathbf{r}).
\end{equation}
The electric field strength $\b E(\b r)$ is complex, the time dependent solution is obtained via $\b E(\b r,t) = \b E(\b r) e^{i\omega t}$. The magnetic field can be obtained by taking a curl of the solution above
\begin{equation}
\mathbf{H}=-\frac{i}{\omega \mu_{0}} \nabla \times \mathbf{E}.
\end{equation}

==== Absorbing boundary conditions ====
Absorbing boundary conditions in the frequency domain can be realized in a number of ways. One of the most straight forward approaches is to impose the Sommerfeld radiation condition or a numerical approximation of it. However, it is more common to surround the computational domain with a thin absorbing layer, which does not reflect outgoing waves. The perfectly matched layer (PML) is one of such absorbing boundary conditions. There are many variants of PML's, but it is clear that the use of one worsens the conditioning of the problem matrix [2]. As it turns out UPML (Uniaxial PML) is numerically less stable than SC-PML (stretched coordinate PML), but efficient preconditioners exist for cases when SC-PML is hard to implement.

===== Uniaxial Perfectly Matched Layer =====
Inside the computational domain we have an amplitude equation for $\b E$,
\begin{equation}
\nabla \times \doubleunderline{\mu}^{-1} \nabla \times \mathbf{E}-\omega^{2} \varepsilon \mathbf{E}=-\mathrm{i} \omega \mathbf{J}.
\end{equation}
We surround the computational domain $\Omega$ with a thin domain $\Omega'$, where we have an uniaxial anisotropic material, which damps the electromagnetic field
\begin{equation}
\nabla \times \doubleunderline{\mu_{s}}^{-1} \nabla \times \mathbf{E}-\omega^{2} \doubleunderline{\varepsilon_{s}} \mathbf{E}=-\mathrm{i} \omega \mathbf{J}
\end{equation}
dielectric and permeability tensors are diagonal,
\begin{equation}
\doubleunderline{\varepsilon_{s}}=\varepsilon\left[\begin{array}{ccc}\frac{s_{y} s_{z}}{s_{x}} & 0 & 0 \\ 0 & \frac{s_{s} s_{x}}{s_{y}} & 0 \\ 0 & 0 & \frac{s_{x} s_{y}}{s_{z}}\end{array}\right], \quad \doubleunderline{\mu_{s}}=\mu\left[\begin{array}{ccc}\frac{s_{y} s_{z}}{s_{x}} & 0 & 0 \\ 0 & \frac{s_{z} s_{x}}{s_{y}} & 0 \\ 0 & 0 & \frac{s_{x} s_{y}}{s_{z}}\end{array}\right].
\end{equation}
To avoid reflections of waves from the boundary between layers, we slowly turn on the loss parameter $s_w$
\begin{equation}
s_{w}(l)=1-\mathrm{i} \frac{\sigma_{w}(l)}{\omega \varepsilon_{0}}
\end{equation}
where $l$ is the distance from the boundary of the PML.

===== Streched Coordinate Perfectly Matched Layer =====
In the stretched coordinate PML the coordinates are stretched and the PML no longer resembles any physical material. The equation inside the PML remains unchanged,
\begin{equation}
\label{eq:SC-PML}
\nabla_{s} \times \mu^{-1} \nabla_{s} \times \mathbf{E}-\omega^{2} \varepsilon \mathbf{E}=-\mathrm{i} \omega \mathbf{J}.
\end{equation}
However, the gradient is redefined as
\begin{equation}
\nabla_{s}=\hat{e}_x \frac{1}{s_{x}} \frac{\partial}{\partial x}+\hat{e}_y \frac{1}{s_{y}} \frac{\partial}{\partial y}+\hat{e}_z \frac{1}{s_{z}} \frac{\partial}{\partial z}.
\end{equation}
Just like with UPML, the damping is slowly turned on. SC-PML is considered more stable than UMPL, but it is usually trickier to implement. For a magnetically homogenous and isotropic medium the SC-PML equation simplifies to
\begin{equation}
-\nabla^2_s\b E + \nabla_s\left( \nabla_{s} \cdot \b E \right) -\omega^2 \eps \mu \b E = -i\omega\b J,
\end{equation}
in matrix form
\begin{multline}
-
\begin{bmatrix}
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_x}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_x}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_x}{z}\right) \\
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_y}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_y}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_y}{z}\right) \\
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_z}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_z}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_z}{z}\right)
\end{bmatrix}
+
\begin{bmatrix}
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_x}{x}\right) +
\frac{1}{s_y s_x}\frac{\partial E_y}{\partial x \partial y} +
\frac{1}{s_z s_x}\frac{\partial E_z}{\partial x \partial z} \\
\frac{1}{s_x s_y}\frac{\partial E_x}{\partial y \partial x} +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_y}{y}\right) +
\frac{1}{s_z s_y}\frac{\partial E_z}{\partial y \partial z} \\
\frac{1}{s_x s_z}\frac{\partial E_x}{\partial z \partial x} +
\frac{1}{s_y s_z}\frac{\partial E_y}{\partial z \partial y} +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_z}{z}\right)
\end{bmatrix} -
\omega^2 \eps \mu \b E
\begin{bmatrix}
E_x \\
E_y \\
E_z
\end{bmatrix}
= -i\omega
\begin{bmatrix}
J_x \\
J_y \\
J_z
\end{bmatrix}
\end{multline}
The expression further simplifies in cases where some of $s_w$'s are constant (cartesian coordinates).

[[File:radial_pml_2D_test_im.png|400px]][[File:radial_pml_2D_test_re.png|400px]][[File:radial_pml_2D_test_sfun.png|400px]]

The above graphs display a simple use case of SC-PML with radial geometry. Instead of being surrounded by a box, the computational domain is surrounded by an annulus of perfectly matched layer, that contains both $s_x$ and $s_y$ components and therefore dampens the wave in both $x$ and $y$ directions throughout the PML. Such approaches are not usually taken in conventional FDTD codes because of difficult implementation, this is not the case when using the Medusa library.

<syntaxhighlight lang="c++" line>
// set equation with source
double l = std::pow(omega/c0, 2.0);
op.lap(source) + l*op.value(source) = -1.0;

for(int i: inter){ // does not contain source node!
op.lap(i) + l*op.value(i) = 0.0;
}
for(int i: PML){ //under assumption sx = sy!
1.0/(sx[i]*sx[i])*op.lap(i) + l*op.value(i) +
((-1.0)/(sx[i]*sx[i]*sx[i])*exop.d1(sx, 0, i)*op.der1(i, 0)
+(-1.0)/(sx[i]*sx[i]*sx[i])*exop.d1(sx, 1, i)*op.der1(i, 1)) = 0.0;
}
for (int i: disc.boundary()){
op.value(i) = 0.0;
}
</syntaxhighlight>

===== Benchmarking =====
Parameters of the PML are usually determined empirically, separately for each simulation, however, there are some rules of thumb. PMLs by Berenger's definition are exact for a continuous equation. When the equation is discretized reflections appear at the boundary of the PML. A lot of care needs to be taken so that the reflections are minimized. This is usually achieved by giving $\sigma_{w}(l)$ a polynomial relationship with the distance $l$ from the PML boundary
\begin{equation}
\sigma_{w}(l,) = \sigma_{max}\left(\frac{l}{d}\right)^m
\end{equation}

where $\sigma_{max}$ is the maximum loss parameter, $d$ is the thickness of the layer and $m$ is the power of the polynomial, and usually lies somewhere between $2$ and $4$. For a normally incident wave on the PML, we can choose the Reflectance of the PML and accordingly tune parameters $d$ and $m$ via
\begin{equation}
\sigma_{w, \max }=-\frac{(m+1) \ln R}{2 \eta_{0} d}
\end{equation}
where $\eta_0$ is the Impedance of free space.

==References==
# Taflove, A. and Hagness, S.C., 2005. Computational electrodynamics: the finite-difference time-domain method. Artech house.
# Shin, W. and Fan, S., 2012. Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell’s equations solvers. Journal of Computational Physics, 231(8), pp.3406-3431.
# Podgornik, Rudolf, and Andrej Vilfan. Elektromagnetno polje. DMFA-založništvo, 2012.
# John, S. G. J., et al. "Photonic crystals: molding the flow of light." Princeton University of press. 2008.

Computational electromagnetics

2020-09-27T08:11:27Z

BlazS: /* Streched Coordinate Perfectly Matched Layer */

== Case studies ==
The following pages describe the basics of computational electromagnetics, starting with relevant derivations and the basics of classical electromagnetism. The subpages include case studies with analytical solutions if they exist and numerical solutions.
* [[Triple dielectric step in 1D]]
* [[Scattering from an infinite cylinder]]
* [[Point source near an anisotropic lens]]

== Classical electromagnetism ==
=== Maxwell's equations in matter ===
Classical electrodynamics is historically one of the most eminent fields of physics as an extension of classical mechanics, it is very successful in explaining a plethora of phenomena. The dynamics of electric and magnetic fields are described with Maxwell's equations. As we will be studying the interaction of electromagnetic waves with different objects, we need Maxwell's equations in matter
<math>
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\begin{align}
&\nabla \times \b{E}(\b{r}, t) = - \dpar{\b{B}(\b{r}, t)}{t}, \label{eq:TFaraday} \\
&\nabla \times \b{H}(\b{r}, t) = \b{j}(\b{r}, t) + \dpar{\b{D}(\b{r}, t)}{t}, \label{eq:TMaxwell-Ampere} \\
&\nabla \cdot \b{D}(\b{r}, t) = \rho(\b{r}, t), \label{eq:TGaussE} \\
&\nabla \cdot \b{B}(\b{r}, t) = 0. \label{eq:TGaussM}
\end{align}
</math>

The system of equations contains four fields. The electric field strength $\b E$ and density $\b D$ and the magnetic field strength $\b H$ and density $\b B$. The four fields are accompanied by the current density $\b j$ and the charge density $\rho$. For a full description of electromagnetic phenomena, we need to provide another two constitutive equations, that relate the strength and density of the fields
<math>
\newcommand{\eps}{\varepsilon}
\begin{align}
\b{B} &= \mu_0 \mu \b{H}, \label{eq:constM} \\
\b{D} &= \eps_0 \eps \b{E} \label{eq:constE}
\end{align}
</math>
where $\varepsilon_0$ and $\mu_0$ are vacuum permittivity and permeability respectively. The dielectric function $\varepsilon$ and magnetic permeability $\mu$ are in general dependant on both $\b E$ and $\b B$ as well as the frequency $\omega$ and can be second order tensors in anisotropic materials. Equations \eqref{eq:constM} and \eqref{eq:constE} already assume linear material properties, generally the polarisation $\b P$ and magnetisation are defined as power series expansions of the electric field density and magnetic field density, as
<math>
\begin{align}
\label{eq:PMexpans}
\b{P} &= \chi_E \b{D} + \mathcal{O}(D^2), \\
\b{M} &= \chi_M \b{H} + \mathcal{O}(H^2).
\end{align}
</math>
The material linearity assumption holds well for small external fields, meaning small $\b D$ and $\b H$. The treatment of nonlinear terms falls within the field of nonlinear optics and is not relevant for our discussion here.

=== Electromagnetic waves ===
<math>
\def\doubleunderline#1{\underline{\underline{#1}}}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\newcommand{\eps}{\varepsilon}
</math>
When studying light and related phenomena, we are usually interested in the behavior of electromagnetic waves, where the electric and magnetic fields oscillate with constant frequency and propagate in a single direction. The wave formulation helps with understanding optical phenomena such as diffraction and reflection and offers easily graspable quantities such as wave amplitude and wavelength. Wave equations for $\b E$ and $\b H$ can be derived from Maxwell's equations for linear and homogeneous materials. Let's start by treating a more general case, where $\doubleunderline \varepsilon$ is a constant second-order tensor and $\mu$ is a scalar quantity. From Maxwell's equations and constitutive relations we obtain
\begin{equation}
\nabla \times \left[ \doubleunderline{\varepsilon}^{-1} \nabla \times \b{H} \right] = \nabla \times \left[ \doubleunderline{\varepsilon}^{-1} \b j \right] - \mu \mu_0 \varepsilon_0 \ddpar{\b H}{t}.
\end{equation}
and
\begin{equation}
\nabla \times \nabla \times \b{E} = \nabla (\nabla \cdot \b E) - \nabla^2 \b E = -\mu\mu_0\dpar{\b j}{t} - \eps_0 \doubleunderline{\varepsilon} \mu \mu_0 \ddpar{\b E}{t}
\end{equation}
which simplify to wave equations for homogenous (linear) materials and no external currents and charges
\begin{equation}
\nabla^2 \b E = \eps\eps_0\mu\mu_0\ddpar{\b E}{t}, \qquad \text{ or } \qquad \nabla^2 \b H = \eps\eps_0\mu\mu_0\ddpar{\b H}{t}.
\end{equation}
One possible way of solving electromagnetic problems is simply to solve the above wave equations numerically. Another approach is to solve Maxwell's equations in the frequency domain. We start by Fourier expansion of the fields $\b E$, $\b D$, $\b B$ and $\b H$ analogous to
\begin{equation}
\label{eq:Fdecomp}
\b{E}(\b{r}, t) = \int\frac{d\omega}{2\pi} \b{E}(\b{r}, \omega) e^{-i\omega t},
\end{equation}
in the frequency domainall the fields become complex. The ratio between the real and complex component of the fields is the phase difference of the material response to external fields. The Fourier expansion leads to the harmonic form of Maxwell's equations
\begin{equation}
\label{eq:FFaraday}
\nabla \times \b{E}(\b{r}, \omega) = - i\omega \b{B} (\b{r}, \omega),
\end{equation}
\begin{equation}
\label{eq:FMaxwell-Ampere}
\nabla \times \b{H}(\b{r}, \omega) = i\omega \b{D}(\b{r}, \omega).
\end{equation}
The constitutive relations are unchanged. If we assume that there are currents we obtain a wave equation for either $\b D$ or $\b H$,
\begin{equation}
\label{eq:frekaniwave}
\nabla \times \left[ \doubleunderline{\eps}^{-1} \nabla \times \b{H} \right] = \omega^2 \mu_0 \eps_0 \mu \b{H}.
\end{equation}
The above form allows for a calculation of a steady state response to an incident harmonic plane wave, meaning that we can study scattering of any incident wave that can be represented by a plane wave expansion. The time dependent solution is obtained by an inverse Fourier transform. In empty space \eqref{eq:frekaniwave} simplifies to the Helmholtz equation
\begin{equation}
\label{eq:frekemptywave}
\nabla^2 \b{H} = - \omega^2 \mu_0 \eps_0 \b{H} = - k^2 \b{H}.
\end{equation}
Wavenumber $k$ in a vacuum is defined as $k = \omega \sqrt{\mu_0 \epsilon_0} = \frac{2 \pi}{\lambda}$. One of the major advantages of the frequency domain approach over the time domain approach is that the frequency dependence of the dielectric function does not need to be considered, as one simply uses the value of the dielectric function at a certain frequency.

===Boundary conditions for Maxwell's equations===
For modeling purposes the treatment of electromagnetic fields on boundaries between materials with different properties, such as $\varepsilon$ and $\mu$, is crucial. The boundary conditions for Maxwell's equations between two linear and homogenous mediums are as follows
\begin{equation}
\begin{array}{rl}
D_{1 n}-D_{2 n}=\sigma, \qquad & E_{1 t}-E_{2 t}=0 \\
B_{1 n}-B_{2 n}=0, \qquad & H_{1 t}-H_{2 t}=K.\end{array}
\end{equation}
Where we have fields $\b{E}_1, \b{D}_1, \b{B}_1, \b{H}_1$ in medium $1$ and fields $\b{E}_2, \b{D}_2, \b{B}_2, \b{H}_2$ in medium $2$. Subscript $t$ and $n$ denote the tangential and normal component respectively, $\sigma$ is the surface charge density and $\b K$ is the surface current density.

===Frequency dependence of the dielectric function===
As mentioned above the dielectric function $\varepsilon$ is a function of frequency. Some analytical properties of the dielectric function $\varepsilon(\omega)$ are well known, due to their relevance in optics, one example being the Kramers-Kronig relations. The dielectric function is in general, based on a nonlocal response function $\chi$, defined as
\begin{equation}
\varepsilon(\omega) = 1 + \int_{0}^{\infty} \chi(\tau) e^{i \omega \tau} \text{d}\tau.
\end{equation}
For a more detailed analysis of the analytical properties of the function, refer to any standard optics or advanced electromagnetics textbook. When it comes to numerical simulations, models of the dielectric function are of more interest than its analytical properties.
====Dielectric materials====
The simplest model for dielectrics is based on the equation of motion for the bound charges
\begin{equation}
m \frac{\mathrm{d}^{2} \mathbf{r}}{\mathrm{d} t^{2}}=-m \gamma \dot{\mathbf{r}}-m \omega_{0}^{2} \mathbf{r}+e \mathbf{E}(t)
\end{equation}
the terms on the right side describe the dissipation of energy, oscillations and the driving force of the external electric field. A short calculation leads to an expression of the frequency dependence of the dielectric function
\begin{equation}
\varepsilon(\omega)=1+ \frac{\omega_p}{\left(\omega_{0}^{2}-\omega^{2}\right)-i \gamma \omega}.
\end{equation}
The above expression can be simplified for high and low frequencies leading to the so-called Debye and Plasmonic relaxation models. The Lorentz model applies for frequencies where no terms in the denominator can be neglected. Typical Debye and Lorentz relaxations are shown in the plots below.

[[File:Debye.png|400px]][[File:Lorentz.png|400px]]

In general, the dielectric function of a material is modeled as a sum of many Debye and Lorentz relaxations, the parameters of which are determined by microscopic properties of the material itself, an extremely interesting case study is the dielectric function of water.
====Conductors====
When dealing with conductors we use the Drude model of conduction
\begin{equation}
m \frac{d \mathbf{v}(t)}{d t}=-m \gamma \mathbf{v}(t)+e \mathbf{E}(t)
\end{equation}
from which we quickly derive the equation for conductivity
\begin{equation}
\sigma(\omega)=\frac{n e^{2}}{m} \frac{1}{\gamma-i \omega},
\end{equation}
which we can relate to the dielectric function
\begin{equation}
\varepsilon(\omega)=1-\frac{\sigma(\omega)}{i \varepsilon_{0} \omega}.
\end{equation}
==Modeling in the time domain==
The single most frequent type of calculation used in computational electromagnetics is the simulation of electromagnetic fields in the time domain. The calculation propagates $\b E(\b r, t)$ and $\b H(\b r, t)$ in time, usually starting with all fields set to zero and some time dependent current source $\b J(\b r, t)$. The main advantage of time domain approaches is the ability to study transient phenomena while its main drawback is the study resonant problems.
===Finite differences in time domain (FDTD)===
==== Formulation ====
Finite differences in time domain were suggested in 1966 by Yee. The method uses finite differences to discretize Maxwell's equations and propagates the fields in time by alternating between updating the electric field strength and magnetic field strength. One of its defining attributes is the use of a staggered grid to represent $\b H$ and $\b E$ nodes. The staggered grid simplifies the calculation of curl expressions and helps with numerical stability and treating the boundaries of materials.

[[File:FDTD2.png|400px|center]]

This type of solution procedure provides an alternative to solving the wave equation for either $\b H$ or $\b E$.

==== Modeling sources ====
Since we are starting with all fields set to zero, we need to somehow introduce power into the model. This is usually done with current sources. The current signal used depends on the type of calculation. If we are interested in a steady-state solution we would use a sinusoidal source, that is slowly turned on and simulated the device until a steady state is reached. On the other hand, if we are interested in calculating the response of a device to many frequencies, we would use a Gaussian pulse and simulate the device until all the energy has left the system.

Taflove in his book [1] differentiates between five types of sources. The easiest to implement are so-called ''hard sources'', where one simply sets the value of a field in a node to the source value e.g.
\begin{equation}
\b E_z(i, t) = \sin 2\pi ft.
\end{equation}
The problem is that the node now acts as a perfect dielectric and incoming waves reflect from the point and do not pass through as they should. On the other hand one can use a ''soft source'' where energy is inserted into the domain via electric and magnetic (non-physical) current
\begin{equation}
\b E_z(i, t) = E_z(i, t) + J(t) \qquad \text{in} \qquad \b H_x(i, t) = H_x(i, t) + M(t).
\end{equation}
Waves no longer reflect from the source, however, the source still radiates in all directions, which is sometimes not desired. This issue is solved with TFSF (total field, scattered field) sources, which are based on the linearity of Maxwell's equations and allow to construct sources that inject the energy only in the desired direction. This is accomplished by decomposing the total fields into the scattered and incident
\begin{equation}
\b{E} = \b{E}^i + \b{E}^s \text{, } \qquad \b{H} = \b{H}^i + \b{H}^s.
\end{equation}
and then dividing the computational domain $\Omega$ into two subdomains. In one of the subdomains, we calculate only using the scattered part of the fields and in the other, we calculate with the total field. With careful treatment of the boundary between the subdomains one can construct an arbitrary incident wave. Alongside the already mentioned source types, ''waveguide'' sources and ''total scattered field'' sources. For an in-depth discussion refer to Taflove [1].

==== Transmittance and Reflectivity ====
Transmittance and reflectivity are dependent on the frequency of the incoming wave. One of the main strengths of time domain calculations is that one can calculate both coefficients for a wide frequency range with a single simulation. This is achieved by using a source that contains a wide range of frequencies i.e. a short Gaussian pulse.

Let us illustrate how to calculate $R$ and $T$ in one dimension. The calculation is similar in higher dimensions. We make two assumptions, firstly the source is directed directly into the device. This means that we will have only the reflected field on one side of the device and only the transmitted field on the other side of the device. Secondly, we assume that nothing reflects from the edges of the domain. We record the signals on both edges of the domain and label them $E_r(t)$ for reflected and $E_t(t)$ for transmitted. It is important to simulate for a long enough time that there is no energy left in the device. Otherwise, the calculation of the coefficients has not converged yet. The Fourier transforms of the signals at the edges, along with the transform of the source, define $R$ and $T$
\begin{equation}
R(\omega) = \left| \frac{E_o(\omega)}{S(\omega)} \right| \qquad \text{and} \qquad T(\omega)= \left| \frac{E_p(\omega)}{S(\omega)} \right|.
\end{equation}
Since we only know the field values at discrete time intervals the transformation is discrete
\begin{equation}
F(\omega) \cong \Delta t \sum_{m=1}^{M}\left(e^{-j \omega \Delta t}\right)^{m} \cdot f(m \Delta t),
\end{equation}
in practice, so that it is not needed to store the signals on both edges which may take up considerable computer memory, the sum can be calculated after each step in the simulation to save space.

<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}
</math>

==Modeling in the frequency domain==
In the frequency domain, two types of calculations are most common. Within the field of photonic crystals, bandstructures of periodic materials are sought by solving a generalized eigenvalue problem $Ax = \omega^2Bx$. In engineering applications, a response of a structure to a periodical source $\b J(x) e^{-i\omega t}$ is more often calculated. We will focus on the second type of calculation.

==== Formulation ====
We need to solve the equation that can be derived by combining Maxwell's equations and constitutive relations
\begin{equation}
\left[(\nabla \times \nabla \times)-\frac{\omega^{2}}{c^{2}} \varepsilon(\mathbf{r})\right] \mathbf{E}(\mathbf{r})=i \omega \mu_{0} \mathbf{J}(\mathbf{r}).
\end{equation}
The electric field strength $\b E(\b r)$ is complex, the time dependent solution is obtained via $\b E(\b r,t) = \b E(\b r) e^{i\omega t}$. The magnetic field can be obtained by taking a curl of the solution above
\begin{equation}
\mathbf{H}=-\frac{i}{\omega \mu_{0}} \nabla \times \mathbf{E}.
\end{equation}

==== Absorbing boundary conditions ====
Absorbing boundary conditions in the frequency domain can be realized in a number of ways. One of the most straight forward approaches is to impose the Sommerfeld radiation condition or a numerical approximation of it. However, it is more common to surround the computational domain with a thin absorbing layer, which does not reflect outgoing waves. The perfectly matched layer (PML) is one of such absorbing boundary conditions. There are many variants of PML's, but it is clear that the use of one worsens the conditioning of the problem matrix [2]. As it turns out UPML (Uniaxial PML) is numerically less stable than SC-PML (stretched coordinate PML), but efficient preconditioners exist for cases when SC-PML is hard to implement.

===== Uniaxial Perfectly Matched Layer =====
Inside the computational domain we have an amplitude equation for $\b E$,
\begin{equation}
\nabla \times \doubleunderline{\mu}^{-1} \nabla \times \mathbf{E}-\omega^{2} \varepsilon \mathbf{E}=-\mathrm{i} \omega \mathbf{J}.
\end{equation}
We surround the computational domain $\Omega$ with a thin domain $\Omega'$, where we have an uniaxial anisotropic material, which damps the electromagnetic field
\begin{equation}
\nabla \times \doubleunderline{\mu_{s}}^{-1} \nabla \times \mathbf{E}-\omega^{2} \doubleunderline{\varepsilon_{s}} \mathbf{E}=-\mathrm{i} \omega \mathbf{J}
\end{equation}
dielectric and permeability tensors are diagonal,
\begin{equation}
\doubleunderline{\varepsilon_{s}}=\varepsilon\left[\begin{array}{ccc}\frac{s_{y} s_{z}}{s_{x}} & 0 & 0 \\ 0 & \frac{s_{s} s_{x}}{s_{y}} & 0 \\ 0 & 0 & \frac{s_{x} s_{y}}{s_{z}}\end{array}\right], \quad \doubleunderline{\mu_{s}}=\mu\left[\begin{array}{ccc}\frac{s_{y} s_{z}}{s_{x}} & 0 & 0 \\ 0 & \frac{s_{z} s_{x}}{s_{y}} & 0 \\ 0 & 0 & \frac{s_{x} s_{y}}{s_{z}}\end{array}\right].
\end{equation}
To avoid reflections of waves from the boundary between layers, we slowly turn on the loss parameter $s_w$
\begin{equation}
s_{w}(l)=1-\mathrm{i} \frac{\sigma_{w}(l)}{\omega \varepsilon_{0}}
\end{equation}
where $l$ is the distance from the boundary of the PML.

===== Streched Coordinate Perfectly Matched Layer =====
In the stretched coordinate PML the coordinates are stretched and the PML no longer resembles any physical material. The equation inside the PML remains unchanged,
\begin{equation}
\label{eq:SC-PML}
\nabla_{s} \times \mu^{-1} \nabla_{s} \times \mathbf{E}-\omega^{2} \varepsilon \mathbf{E}=-\mathrm{i} \omega \mathbf{J}.
\end{equation}
However, the gradient is redefined as
\begin{equation}
\nabla_{s}=\hat{e}_x \frac{1}{s_{x}} \frac{\partial}{\partial x}+\hat{e}_y \frac{1}{s_{y}} \frac{\partial}{\partial y}+\hat{e}_z \frac{1}{s_{z}} \frac{\partial}{\partial z}.
\end{equation}
Just like with UPML, the damping is slowly turned on. SC-PML is considered more stable than UMPL, but it is usually trickier to implement. For a magnetically homogenous and isotropic medium the SC-PML equation simplifies to
\begin{equation}
-\nabla^2_s\b E + \nabla_s\left( \nabla_{s} \cdot \b E \right) -\omega^2 \eps \mu \b E = -i\omega\b J,
\end{equation}
in matrix form
\begin{multline}
-
\begin{bmatrix}
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_x}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_x}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_x}{z}\right) \\
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_y}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_y}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_y}{z}\right) \\
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_z}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_z}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_z}{z}\right)
\end{bmatrix}
+
\begin{bmatrix}
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_x}{x}\right) +
\frac{1}{s_y s_x}\frac{\partial E_y}{\partial x \partial y} +
\frac{1}{s_z s_x}\frac{\partial E_z}{\partial z \partial z} \\
\frac{1}{s_x s_y}\frac{\partial E_x}{\partial y \partial x} +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_y}{y}\right) +
\frac{1}{s_z s_y}\frac{\partial E_z}{\partial y \partial z} \\
\frac{1}{s_x s_z}\frac{\partial E_x}{\partial z \partial x} +
\frac{1}{s_y s_z}\frac{\partial E_y}{\partial z \partial y} +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_z}{z}\right)
\end{bmatrix} -
\omega^2 \eps \mu \b E
\begin{bmatrix}
E_x \\
E_y \\
E_z
\end{bmatrix}
= -i\omega
\begin{bmatrix}
J_x \\
J_y \\
J_z
\end{bmatrix}
\end{multline}
The expression further simplifies in cases where some of $s_w$'s are constant (cartesian coordinates).

[[File:radial_pml_2D_test_im.png|400px]][[File:radial_pml_2D_test_re.png|400px]][[File:radial_pml_2D_test_sfun.png|400px]]

The above graphs display a simple use case of SC-PML with radial geometry. Instead of being surrounded by a box, the computational domain is surrounded by an annulus of perfectly matched layer, that contains both $s_x$ and $s_y$ components and therefore dampens the wave in both $x$ and $y$ directions throughout the PML. Such approaches are not usually taken in conventional FDTD codes because of difficult implementation, this is not the case when using the Medusa library.

<syntaxhighlight lang="c++" line>
// set equation with source
double l = std::pow(omega/c0, 2.0);
op.lap(source) + l*op.value(source) = -1.0;

for(int i: inter){ // does not contain source node!
op.lap(i) + l*op.value(i) = 0.0;
}
for(int i: PML){ //under assumption sx = sy!
1.0/(sx[i]*sx[i])*op.lap(i) + l*op.value(i) +
((-1.0)/(sx[i]*sx[i]*sx[i])*exop.d1(sx, 0, i)*op.der1(i, 0)
+(-1.0)/(sx[i]*sx[i]*sx[i])*exop.d1(sx, 1, i)*op.der1(i, 1)) = 0.0;
}
for (int i: disc.boundary()){
op.value(i) = 0.0;
}
</syntaxhighlight>

===== Benchmarking =====
Parameters of the PML are usually determined empirically, separately for each simulation, however, there are some rules of thumb. PMLs by Berenger's definition are exact for a continuous equation. When the equation is discretized reflections appear at the boundary of the PML. A lot of care needs to be taken so that the reflections are minimized. This is usually achieved by giving $\sigma_{w}(l)$ a polynomial relationship with the distance $l$ from the PML boundary
\begin{equation}
\sigma_{w}(l,) = \sigma_{max}\left(\frac{l}{d}\right)^m
\end{equation}

where $\sigma_{max}$ is the maximum loss parameter, $d$ is the thickness of the layer and $m$ is the power of the polynomial, and usually lies somewhere between $2$ and $4$. For a normally incident wave on the PML, we can choose the Reflectance of the PML and accordingly tune parameters $d$ and $m$ via
\begin{equation}
\sigma_{w, \max }=-\frac{(m+1) \ln R}{2 \eta_{0} d}
\end{equation}
where $\eta_0$ is the Impedance of free space.

==References==
# Taflove, A. and Hagness, S.C., 2005. Computational electrodynamics: the finite-difference time-domain method. Artech house.
# Shin, W. and Fan, S., 2012. Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell’s equations solvers. Journal of Computational Physics, 231(8), pp.3406-3431.
# Podgornik, Rudolf, and Andrej Vilfan. Elektromagnetno polje. DMFA-založništvo, 2012.
# John, S. G. J., et al. "Photonic crystals: molding the flow of light." Princeton University of press. 2008.

Computational electromagnetics

2020-09-27T08:10:14Z

BlazS: /* Streched Coordinate Perfectly Matched Layer */

== Case studies ==
The following pages describe the basics of computational electromagnetics, starting with relevant derivations and the basics of classical electromagnetism. The subpages include case studies with analytical solutions if they exist and numerical solutions.
* [[Triple dielectric step in 1D]]
* [[Scattering from an infinite cylinder]]
* [[Point source near an anisotropic lens]]

== Classical electromagnetism ==
=== Maxwell's equations in matter ===
Classical electrodynamics is historically one of the most eminent fields of physics as an extension of classical mechanics, it is very successful in explaining a plethora of phenomena. The dynamics of electric and magnetic fields are described with Maxwell's equations. As we will be studying the interaction of electromagnetic waves with different objects, we need Maxwell's equations in matter
<math>
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\begin{align}
&\nabla \times \b{E}(\b{r}, t) = - \dpar{\b{B}(\b{r}, t)}{t}, \label{eq:TFaraday} \\
&\nabla \times \b{H}(\b{r}, t) = \b{j}(\b{r}, t) + \dpar{\b{D}(\b{r}, t)}{t}, \label{eq:TMaxwell-Ampere} \\
&\nabla \cdot \b{D}(\b{r}, t) = \rho(\b{r}, t), \label{eq:TGaussE} \\
&\nabla \cdot \b{B}(\b{r}, t) = 0. \label{eq:TGaussM}
\end{align}
</math>

The system of equations contains four fields. The electric field strength $\b E$ and density $\b D$ and the magnetic field strength $\b H$ and density $\b B$. The four fields are accompanied by the current density $\b j$ and the charge density $\rho$. For a full description of electromagnetic phenomena, we need to provide another two constitutive equations, that relate the strength and density of the fields
<math>
\newcommand{\eps}{\varepsilon}
\begin{align}
\b{B} &= \mu_0 \mu \b{H}, \label{eq:constM} \\
\b{D} &= \eps_0 \eps \b{E} \label{eq:constE}
\end{align}
</math>
where $\varepsilon_0$ and $\mu_0$ are vacuum permittivity and permeability respectively. The dielectric function $\varepsilon$ and magnetic permeability $\mu$ are in general dependant on both $\b E$ and $\b B$ as well as the frequency $\omega$ and can be second order tensors in anisotropic materials. Equations \eqref{eq:constM} and \eqref{eq:constE} already assume linear material properties, generally the polarisation $\b P$ and magnetisation are defined as power series expansions of the electric field density and magnetic field density, as
<math>
\begin{align}
\label{eq:PMexpans}
\b{P} &= \chi_E \b{D} + \mathcal{O}(D^2), \\
\b{M} &= \chi_M \b{H} + \mathcal{O}(H^2).
\end{align}
</math>
The material linearity assumption holds well for small external fields, meaning small $\b D$ and $\b H$. The treatment of nonlinear terms falls within the field of nonlinear optics and is not relevant for our discussion here.

=== Electromagnetic waves ===
<math>
\def\doubleunderline#1{\underline{\underline{#1}}}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\newcommand{\eps}{\varepsilon}
</math>
When studying light and related phenomena, we are usually interested in the behavior of electromagnetic waves, where the electric and magnetic fields oscillate with constant frequency and propagate in a single direction. The wave formulation helps with understanding optical phenomena such as diffraction and reflection and offers easily graspable quantities such as wave amplitude and wavelength. Wave equations for $\b E$ and $\b H$ can be derived from Maxwell's equations for linear and homogeneous materials. Let's start by treating a more general case, where $\doubleunderline \varepsilon$ is a constant second-order tensor and $\mu$ is a scalar quantity. From Maxwell's equations and constitutive relations we obtain
\begin{equation}
\nabla \times \left[ \doubleunderline{\varepsilon}^{-1} \nabla \times \b{H} \right] = \nabla \times \left[ \doubleunderline{\varepsilon}^{-1} \b j \right] - \mu \mu_0 \varepsilon_0 \ddpar{\b H}{t}.
\end{equation}
and
\begin{equation}
\nabla \times \nabla \times \b{E} = \nabla (\nabla \cdot \b E) - \nabla^2 \b E = -\mu\mu_0\dpar{\b j}{t} - \eps_0 \doubleunderline{\varepsilon} \mu \mu_0 \ddpar{\b E}{t}
\end{equation}
which simplify to wave equations for homogenous (linear) materials and no external currents and charges
\begin{equation}
\nabla^2 \b E = \eps\eps_0\mu\mu_0\ddpar{\b E}{t}, \qquad \text{ or } \qquad \nabla^2 \b H = \eps\eps_0\mu\mu_0\ddpar{\b H}{t}.
\end{equation}
One possible way of solving electromagnetic problems is simply to solve the above wave equations numerically. Another approach is to solve Maxwell's equations in the frequency domain. We start by Fourier expansion of the fields $\b E$, $\b D$, $\b B$ and $\b H$ analogous to
\begin{equation}
\label{eq:Fdecomp}
\b{E}(\b{r}, t) = \int\frac{d\omega}{2\pi} \b{E}(\b{r}, \omega) e^{-i\omega t},
\end{equation}
in the frequency domainall the fields become complex. The ratio between the real and complex component of the fields is the phase difference of the material response to external fields. The Fourier expansion leads to the harmonic form of Maxwell's equations
\begin{equation}
\label{eq:FFaraday}
\nabla \times \b{E}(\b{r}, \omega) = - i\omega \b{B} (\b{r}, \omega),
\end{equation}
\begin{equation}
\label{eq:FMaxwell-Ampere}
\nabla \times \b{H}(\b{r}, \omega) = i\omega \b{D}(\b{r}, \omega).
\end{equation}
The constitutive relations are unchanged. If we assume that there are currents we obtain a wave equation for either $\b D$ or $\b H$,
\begin{equation}
\label{eq:frekaniwave}
\nabla \times \left[ \doubleunderline{\eps}^{-1} \nabla \times \b{H} \right] = \omega^2 \mu_0 \eps_0 \mu \b{H}.
\end{equation}
The above form allows for a calculation of a steady state response to an incident harmonic plane wave, meaning that we can study scattering of any incident wave that can be represented by a plane wave expansion. The time dependent solution is obtained by an inverse Fourier transform. In empty space \eqref{eq:frekaniwave} simplifies to the Helmholtz equation
\begin{equation}
\label{eq:frekemptywave}
\nabla^2 \b{H} = - \omega^2 \mu_0 \eps_0 \b{H} = - k^2 \b{H}.
\end{equation}
Wavenumber $k$ in a vacuum is defined as $k = \omega \sqrt{\mu_0 \epsilon_0} = \frac{2 \pi}{\lambda}$. One of the major advantages of the frequency domain approach over the time domain approach is that the frequency dependence of the dielectric function does not need to be considered, as one simply uses the value of the dielectric function at a certain frequency.

===Boundary conditions for Maxwell's equations===
For modeling purposes the treatment of electromagnetic fields on boundaries between materials with different properties, such as $\varepsilon$ and $\mu$, is crucial. The boundary conditions for Maxwell's equations between two linear and homogenous mediums are as follows
\begin{equation}
\begin{array}{rl}
D_{1 n}-D_{2 n}=\sigma, \qquad & E_{1 t}-E_{2 t}=0 \\
B_{1 n}-B_{2 n}=0, \qquad & H_{1 t}-H_{2 t}=K.\end{array}
\end{equation}
Where we have fields $\b{E}_1, \b{D}_1, \b{B}_1, \b{H}_1$ in medium $1$ and fields $\b{E}_2, \b{D}_2, \b{B}_2, \b{H}_2$ in medium $2$. Subscript $t$ and $n$ denote the tangential and normal component respectively, $\sigma$ is the surface charge density and $\b K$ is the surface current density.

===Frequency dependence of the dielectric function===
As mentioned above the dielectric function $\varepsilon$ is a function of frequency. Some analytical properties of the dielectric function $\varepsilon(\omega)$ are well known, due to their relevance in optics, one example being the Kramers-Kronig relations. The dielectric function is in general, based on a nonlocal response function $\chi$, defined as
\begin{equation}
\varepsilon(\omega) = 1 + \int_{0}^{\infty} \chi(\tau) e^{i \omega \tau} \text{d}\tau.
\end{equation}
For a more detailed analysis of the analytical properties of the function, refer to any standard optics or advanced electromagnetics textbook. When it comes to numerical simulations, models of the dielectric function are of more interest than its analytical properties.
====Dielectric materials====
The simplest model for dielectrics is based on the equation of motion for the bound charges
\begin{equation}
m \frac{\mathrm{d}^{2} \mathbf{r}}{\mathrm{d} t^{2}}=-m \gamma \dot{\mathbf{r}}-m \omega_{0}^{2} \mathbf{r}+e \mathbf{E}(t)
\end{equation}
the terms on the right side describe the dissipation of energy, oscillations and the driving force of the external electric field. A short calculation leads to an expression of the frequency dependence of the dielectric function
\begin{equation}
\varepsilon(\omega)=1+ \frac{\omega_p}{\left(\omega_{0}^{2}-\omega^{2}\right)-i \gamma \omega}.
\end{equation}
The above expression can be simplified for high and low frequencies leading to the so-called Debye and Plasmonic relaxation models. The Lorentz model applies for frequencies where no terms in the denominator can be neglected. Typical Debye and Lorentz relaxations are shown in the plots below.

[[File:Debye.png|400px]][[File:Lorentz.png|400px]]

In general, the dielectric function of a material is modeled as a sum of many Debye and Lorentz relaxations, the parameters of which are determined by microscopic properties of the material itself, an extremely interesting case study is the dielectric function of water.
====Conductors====
When dealing with conductors we use the Drude model of conduction
\begin{equation}
m \frac{d \mathbf{v}(t)}{d t}=-m \gamma \mathbf{v}(t)+e \mathbf{E}(t)
\end{equation}
from which we quickly derive the equation for conductivity
\begin{equation}
\sigma(\omega)=\frac{n e^{2}}{m} \frac{1}{\gamma-i \omega},
\end{equation}
which we can relate to the dielectric function
\begin{equation}
\varepsilon(\omega)=1-\frac{\sigma(\omega)}{i \varepsilon_{0} \omega}.
\end{equation}
==Modeling in the time domain==
The single most frequent type of calculation used in computational electromagnetics is the simulation of electromagnetic fields in the time domain. The calculation propagates $\b E(\b r, t)$ and $\b H(\b r, t)$ in time, usually starting with all fields set to zero and some time dependent current source $\b J(\b r, t)$. The main advantage of time domain approaches is the ability to study transient phenomena while its main drawback is the study resonant problems.
===Finite differences in time domain (FDTD)===
==== Formulation ====
Finite differences in time domain were suggested in 1966 by Yee. The method uses finite differences to discretize Maxwell's equations and propagates the fields in time by alternating between updating the electric field strength and magnetic field strength. One of its defining attributes is the use of a staggered grid to represent $\b H$ and $\b E$ nodes. The staggered grid simplifies the calculation of curl expressions and helps with numerical stability and treating the boundaries of materials.

[[File:FDTD2.png|400px|center]]

This type of solution procedure provides an alternative to solving the wave equation for either $\b H$ or $\b E$.

==== Modeling sources ====
Since we are starting with all fields set to zero, we need to somehow introduce power into the model. This is usually done with current sources. The current signal used depends on the type of calculation. If we are interested in a steady-state solution we would use a sinusoidal source, that is slowly turned on and simulated the device until a steady state is reached. On the other hand, if we are interested in calculating the response of a device to many frequencies, we would use a Gaussian pulse and simulate the device until all the energy has left the system.

Taflove in his book [1] differentiates between five types of sources. The easiest to implement are so-called ''hard sources'', where one simply sets the value of a field in a node to the source value e.g.
\begin{equation}
\b E_z(i, t) = \sin 2\pi ft.
\end{equation}
The problem is that the node now acts as a perfect dielectric and incoming waves reflect from the point and do not pass through as they should. On the other hand one can use a ''soft source'' where energy is inserted into the domain via electric and magnetic (non-physical) current
\begin{equation}
\b E_z(i, t) = E_z(i, t) + J(t) \qquad \text{in} \qquad \b H_x(i, t) = H_x(i, t) + M(t).
\end{equation}
Waves no longer reflect from the source, however, the source still radiates in all directions, which is sometimes not desired. This issue is solved with TFSF (total field, scattered field) sources, which are based on the linearity of Maxwell's equations and allow to construct sources that inject the energy only in the desired direction. This is accomplished by decomposing the total fields into the scattered and incident
\begin{equation}
\b{E} = \b{E}^i + \b{E}^s \text{, } \qquad \b{H} = \b{H}^i + \b{H}^s.
\end{equation}
and then dividing the computational domain $\Omega$ into two subdomains. In one of the subdomains, we calculate only using the scattered part of the fields and in the other, we calculate with the total field. With careful treatment of the boundary between the subdomains one can construct an arbitrary incident wave. Alongside the already mentioned source types, ''waveguide'' sources and ''total scattered field'' sources. For an in-depth discussion refer to Taflove [1].

==== Transmittance and Reflectivity ====
Transmittance and reflectivity are dependent on the frequency of the incoming wave. One of the main strengths of time domain calculations is that one can calculate both coefficients for a wide frequency range with a single simulation. This is achieved by using a source that contains a wide range of frequencies i.e. a short Gaussian pulse.

Let us illustrate how to calculate $R$ and $T$ in one dimension. The calculation is similar in higher dimensions. We make two assumptions, firstly the source is directed directly into the device. This means that we will have only the reflected field on one side of the device and only the transmitted field on the other side of the device. Secondly, we assume that nothing reflects from the edges of the domain. We record the signals on both edges of the domain and label them $E_r(t)$ for reflected and $E_t(t)$ for transmitted. It is important to simulate for a long enough time that there is no energy left in the device. Otherwise, the calculation of the coefficients has not converged yet. The Fourier transforms of the signals at the edges, along with the transform of the source, define $R$ and $T$
\begin{equation}
R(\omega) = \left| \frac{E_o(\omega)}{S(\omega)} \right| \qquad \text{and} \qquad T(\omega)= \left| \frac{E_p(\omega)}{S(\omega)} \right|.
\end{equation}
Since we only know the field values at discrete time intervals the transformation is discrete
\begin{equation}
F(\omega) \cong \Delta t \sum_{m=1}^{M}\left(e^{-j \omega \Delta t}\right)^{m} \cdot f(m \Delta t),
\end{equation}
in practice, so that it is not needed to store the signals on both edges which may take up considerable computer memory, the sum can be calculated after each step in the simulation to save space.

<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}
</math>

==Modeling in the frequency domain==
In the frequency domain, two types of calculations are most common. Within the field of photonic crystals, bandstructures of periodic materials are sought by solving a generalized eigenvalue problem $Ax = \omega^2Bx$. In engineering applications, a response of a structure to a periodical source $\b J(x) e^{-i\omega t}$ is more often calculated. We will focus on the second type of calculation.

==== Formulation ====
We need to solve the equation that can be derived by combining Maxwell's equations and constitutive relations
\begin{equation}
\left[(\nabla \times \nabla \times)-\frac{\omega^{2}}{c^{2}} \varepsilon(\mathbf{r})\right] \mathbf{E}(\mathbf{r})=i \omega \mu_{0} \mathbf{J}(\mathbf{r}).
\end{equation}
The electric field strength $\b E(\b r)$ is complex, the time dependent solution is obtained via $\b E(\b r,t) = \b E(\b r) e^{i\omega t}$. The magnetic field can be obtained by taking a curl of the solution above
\begin{equation}
\mathbf{H}=-\frac{i}{\omega \mu_{0}} \nabla \times \mathbf{E}.
\end{equation}

==== Absorbing boundary conditions ====
Absorbing boundary conditions in the frequency domain can be realized in a number of ways. One of the most straight forward approaches is to impose the Sommerfeld radiation condition or a numerical approximation of it. However, it is more common to surround the computational domain with a thin absorbing layer, which does not reflect outgoing waves. The perfectly matched layer (PML) is one of such absorbing boundary conditions. There are many variants of PML's, but it is clear that the use of one worsens the conditioning of the problem matrix [2]. As it turns out UPML (Uniaxial PML) is numerically less stable than SC-PML (stretched coordinate PML), but efficient preconditioners exist for cases when SC-PML is hard to implement.

===== Uniaxial Perfectly Matched Layer =====
Inside the computational domain we have an amplitude equation for $\b E$,
\begin{equation}
\nabla \times \doubleunderline{\mu}^{-1} \nabla \times \mathbf{E}-\omega^{2} \varepsilon \mathbf{E}=-\mathrm{i} \omega \mathbf{J}.
\end{equation}
We surround the computational domain $\Omega$ with a thin domain $\Omega'$, where we have an uniaxial anisotropic material, which damps the electromagnetic field
\begin{equation}
\nabla \times \doubleunderline{\mu_{s}}^{-1} \nabla \times \mathbf{E}-\omega^{2} \doubleunderline{\varepsilon_{s}} \mathbf{E}=-\mathrm{i} \omega \mathbf{J}
\end{equation}
dielectric and permeability tensors are diagonal,
\begin{equation}
\doubleunderline{\varepsilon_{s}}=\varepsilon\left[\begin{array}{ccc}\frac{s_{y} s_{z}}{s_{x}} & 0 & 0 \\ 0 & \frac{s_{s} s_{x}}{s_{y}} & 0 \\ 0 & 0 & \frac{s_{x} s_{y}}{s_{z}}\end{array}\right], \quad \doubleunderline{\mu_{s}}=\mu\left[\begin{array}{ccc}\frac{s_{y} s_{z}}{s_{x}} & 0 & 0 \\ 0 & \frac{s_{z} s_{x}}{s_{y}} & 0 \\ 0 & 0 & \frac{s_{x} s_{y}}{s_{z}}\end{array}\right].
\end{equation}
To avoid reflections of waves from the boundary between layers, we slowly turn on the loss parameter $s_w$
\begin{equation}
s_{w}(l)=1-\mathrm{i} \frac{\sigma_{w}(l)}{\omega \varepsilon_{0}}
\end{equation}
where $l$ is the distance from the boundary of the PML.

===== Streched Coordinate Perfectly Matched Layer =====
In the stretched coordinate PML the coordinates are stretched and the PML no longer resembles any physical material. The equation inside the PML remains unchanged,
\begin{equation}
\label{eq:SC-PML}
\nabla_{s} \times \mu^{-1} \nabla_{s} \times \mathbf{E}-\omega^{2} \varepsilon \mathbf{E}=-\mathrm{i} \omega \mathbf{J}.
\end{equation}
However, the gradient is redefined as
\begin{equation}
\nabla_{s}=\hat{e}_x \frac{1}{s_{x}} \frac{\partial}{\partial x}+\hat{e}_y \frac{1}{s_{y}} \frac{\partial}{\partial y}+\hat{e}_z \frac{1}{s_{z}} \frac{\partial}{\partial z}.
\end{equation}
Just like with UPML, the damping is slowly turned on. SC-PML is considered more stable than UMPL, but it is usually trickier to implement. For a magnetically homogenous and isotropic medium the SC-PML equation simplifies to
\begin{equation}
-\nabla^2_s\b E + \nabla_s\left( \nabla_{s} \cdot \b E \right) -\omega^2 \eps \mu \b E = -i\omega\b J,
\end{equation}
in matrix form
\begin{multline}
-
\begin{bmatrix}
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_x}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_x}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_x}{z}\right) \\
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_y}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_y}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_y}{z}\right) \\
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_z}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_z}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_z}{z}\right)
\end{bmatrix}
+
\begin{bmatrix}
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_x}{x}\right) +
\frac{1}{s_y s_x}\frac{\partial E_y}{\partial x \partial y} +
\frac{1}{s_z s_x}\frac{\partial E_z}{\partial z \partial z} \\
\frac{1}{s_x s_y}\frac{\partial E_x}{\partial y \partial x} +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_y}{y}\right) +
\frac{1}{s_z s_y}\frac{\partial E_z}{\partial y \partial z} \\
\frac{1}{s_x s_z}\frac{\partial E_x}{\partial z \partial x} +
\frac{1}{s_y s_z}\frac{\partial E_y}{\partial z \partial y} +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_z}{z}\right)
\end{bmatrix} -
\\ -\omega^2 \eps \mu \b E
\begin{bmatrix}
E_x \\
E_y \\
E_z
\end{bmatrix}
= -i\omega
\begin{bmatrix}
J_x \\
J_y \\
J_z
\end{bmatrix}
\end{multline}
The expression further simplifies in cases where some of $s_w$'s are constant (cartesian coordinates).

[[File:radial_pml_2D_test_im.png|400px]][[File:radial_pml_2D_test_re.png|400px]][[File:radial_pml_2D_test_sfun.png|400px]]

The above graphs display a simple use case of SC-PML with radial geometry. Instead of being surrounded by a box, the computational domain is surrounded by an annulus of perfectly matched layer, that contains both $s_x$ and $s_y$ components and therefore dampens the wave in both $x$ and $y$ directions throughout the PML. Such approaches are not usually taken in conventional FDTD codes because of difficult implementation, this is not the case when using the Medusa library.

<syntaxhighlight lang="c++" line>
// set equation with source
double l = std::pow(omega/c0, 2.0);
op.lap(source) + l*op.value(source) = -1.0;

for(int i: inter){ // does not contain source node!
op.lap(i) + l*op.value(i) = 0.0;
}
for(int i: PML){ //under assumption sx = sy!
1.0/(sx[i]*sx[i])*op.lap(i) + l*op.value(i) +
((-1.0)/(sx[i]*sx[i]*sx[i])*exop.d1(sx, 0, i)*op.der1(i, 0)
+(-1.0)/(sx[i]*sx[i]*sx[i])*exop.d1(sx, 1, i)*op.der1(i, 1)) = 0.0;
}
for (int i: disc.boundary()){
op.value(i) = 0.0;
}
</syntaxhighlight>

===== Benchmarking =====
Parameters of the PML are usually determined empirically, separately for each simulation, however, there are some rules of thumb. PMLs by Berenger's definition are exact for a continuous equation. When the equation is discretized reflections appear at the boundary of the PML. A lot of care needs to be taken so that the reflections are minimized. This is usually achieved by giving $\sigma_{w}(l)$ a polynomial relationship with the distance $l$ from the PML boundary
\begin{equation}
\sigma_{w}(l,) = \sigma_{max}\left(\frac{l}{d}\right)^m
\end{equation}

where $\sigma_{max}$ is the maximum loss parameter, $d$ is the thickness of the layer and $m$ is the power of the polynomial, and usually lies somewhere between $2$ and $4$. For a normally incident wave on the PML, we can choose the Reflectance of the PML and accordingly tune parameters $d$ and $m$ via
\begin{equation}
\sigma_{w, \max }=-\frac{(m+1) \ln R}{2 \eta_{0} d}
\end{equation}
where $\eta_0$ is the Impedance of free space.

==References==
# Taflove, A. and Hagness, S.C., 2005. Computational electrodynamics: the finite-difference time-domain method. Artech house.
# Shin, W. and Fan, S., 2012. Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell’s equations solvers. Journal of Computational Physics, 231(8), pp.3406-3431.
# Podgornik, Rudolf, and Andrej Vilfan. Elektromagnetno polje. DMFA-založništvo, 2012.
# John, S. G. J., et al. "Photonic crystals: molding the flow of light." Princeton University of press. 2008.

Computational electromagnetics

2020-09-27T08:07:27Z

BlazS: /* Streched Coordinate Perfectly Matched Layer */

== Case studies ==
The following pages describe the basics of computational electromagnetics, starting with relevant derivations and the basics of classical electromagnetism. The subpages include case studies with analytical solutions if they exist and numerical solutions.
* [[Triple dielectric step in 1D]]
* [[Scattering from an infinite cylinder]]
* [[Point source near an anisotropic lens]]

== Classical electromagnetism ==
=== Maxwell's equations in matter ===
Classical electrodynamics is historically one of the most eminent fields of physics as an extension of classical mechanics, it is very successful in explaining a plethora of phenomena. The dynamics of electric and magnetic fields are described with Maxwell's equations. As we will be studying the interaction of electromagnetic waves with different objects, we need Maxwell's equations in matter
<math>
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\begin{align}
&\nabla \times \b{E}(\b{r}, t) = - \dpar{\b{B}(\b{r}, t)}{t}, \label{eq:TFaraday} \\
&\nabla \times \b{H}(\b{r}, t) = \b{j}(\b{r}, t) + \dpar{\b{D}(\b{r}, t)}{t}, \label{eq:TMaxwell-Ampere} \\
&\nabla \cdot \b{D}(\b{r}, t) = \rho(\b{r}, t), \label{eq:TGaussE} \\
&\nabla \cdot \b{B}(\b{r}, t) = 0. \label{eq:TGaussM}
\end{align}
</math>

The system of equations contains four fields. The electric field strength $\b E$ and density $\b D$ and the magnetic field strength $\b H$ and density $\b B$. The four fields are accompanied by the current density $\b j$ and the charge density $\rho$. For a full description of electromagnetic phenomena, we need to provide another two constitutive equations, that relate the strength and density of the fields
<math>
\newcommand{\eps}{\varepsilon}
\begin{align}
\b{B} &= \mu_0 \mu \b{H}, \label{eq:constM} \\
\b{D} &= \eps_0 \eps \b{E} \label{eq:constE}
\end{align}
</math>
where $\varepsilon_0$ and $\mu_0$ are vacuum permittivity and permeability respectively. The dielectric function $\varepsilon$ and magnetic permeability $\mu$ are in general dependant on both $\b E$ and $\b B$ as well as the frequency $\omega$ and can be second order tensors in anisotropic materials. Equations \eqref{eq:constM} and \eqref{eq:constE} already assume linear material properties, generally the polarisation $\b P$ and magnetisation are defined as power series expansions of the electric field density and magnetic field density, as
<math>
\begin{align}
\label{eq:PMexpans}
\b{P} &= \chi_E \b{D} + \mathcal{O}(D^2), \\
\b{M} &= \chi_M \b{H} + \mathcal{O}(H^2).
\end{align}
</math>
The material linearity assumption holds well for small external fields, meaning small $\b D$ and $\b H$. The treatment of nonlinear terms falls within the field of nonlinear optics and is not relevant for our discussion here.

=== Electromagnetic waves ===
<math>
\def\doubleunderline#1{\underline{\underline{#1}}}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\newcommand{\eps}{\varepsilon}
</math>
When studying light and related phenomena, we are usually interested in the behavior of electromagnetic waves, where the electric and magnetic fields oscillate with constant frequency and propagate in a single direction. The wave formulation helps with understanding optical phenomena such as diffraction and reflection and offers easily graspable quantities such as wave amplitude and wavelength. Wave equations for $\b E$ and $\b H$ can be derived from Maxwell's equations for linear and homogeneous materials. Let's start by treating a more general case, where $\doubleunderline \varepsilon$ is a constant second-order tensor and $\mu$ is a scalar quantity. From Maxwell's equations and constitutive relations we obtain
\begin{equation}
\nabla \times \left[ \doubleunderline{\varepsilon}^{-1} \nabla \times \b{H} \right] = \nabla \times \left[ \doubleunderline{\varepsilon}^{-1} \b j \right] - \mu \mu_0 \varepsilon_0 \ddpar{\b H}{t}.
\end{equation}
and
\begin{equation}
\nabla \times \nabla \times \b{E} = \nabla (\nabla \cdot \b E) - \nabla^2 \b E = -\mu\mu_0\dpar{\b j}{t} - \eps_0 \doubleunderline{\varepsilon} \mu \mu_0 \ddpar{\b E}{t}
\end{equation}
which simplify to wave equations for homogenous (linear) materials and no external currents and charges
\begin{equation}
\nabla^2 \b E = \eps\eps_0\mu\mu_0\ddpar{\b E}{t}, \qquad \text{ or } \qquad \nabla^2 \b H = \eps\eps_0\mu\mu_0\ddpar{\b H}{t}.
\end{equation}
One possible way of solving electromagnetic problems is simply to solve the above wave equations numerically. Another approach is to solve Maxwell's equations in the frequency domain. We start by Fourier expansion of the fields $\b E$, $\b D$, $\b B$ and $\b H$ analogous to
\begin{equation}
\label{eq:Fdecomp}
\b{E}(\b{r}, t) = \int\frac{d\omega}{2\pi} \b{E}(\b{r}, \omega) e^{-i\omega t},
\end{equation}
in the frequency domainall the fields become complex. The ratio between the real and complex component of the fields is the phase difference of the material response to external fields. The Fourier expansion leads to the harmonic form of Maxwell's equations
\begin{equation}
\label{eq:FFaraday}
\nabla \times \b{E}(\b{r}, \omega) = - i\omega \b{B} (\b{r}, \omega),
\end{equation}
\begin{equation}
\label{eq:FMaxwell-Ampere}
\nabla \times \b{H}(\b{r}, \omega) = i\omega \b{D}(\b{r}, \omega).
\end{equation}
The constitutive relations are unchanged. If we assume that there are currents we obtain a wave equation for either $\b D$ or $\b H$,
\begin{equation}
\label{eq:frekaniwave}
\nabla \times \left[ \doubleunderline{\eps}^{-1} \nabla \times \b{H} \right] = \omega^2 \mu_0 \eps_0 \mu \b{H}.
\end{equation}
The above form allows for a calculation of a steady state response to an incident harmonic plane wave, meaning that we can study scattering of any incident wave that can be represented by a plane wave expansion. The time dependent solution is obtained by an inverse Fourier transform. In empty space \eqref{eq:frekaniwave} simplifies to the Helmholtz equation
\begin{equation}
\label{eq:frekemptywave}
\nabla^2 \b{H} = - \omega^2 \mu_0 \eps_0 \b{H} = - k^2 \b{H}.
\end{equation}
Wavenumber $k$ in a vacuum is defined as $k = \omega \sqrt{\mu_0 \epsilon_0} = \frac{2 \pi}{\lambda}$. One of the major advantages of the frequency domain approach over the time domain approach is that the frequency dependence of the dielectric function does not need to be considered, as one simply uses the value of the dielectric function at a certain frequency.

===Boundary conditions for Maxwell's equations===
For modeling purposes the treatment of electromagnetic fields on boundaries between materials with different properties, such as $\varepsilon$ and $\mu$, is crucial. The boundary conditions for Maxwell's equations between two linear and homogenous mediums are as follows
\begin{equation}
\begin{array}{rl}
D_{1 n}-D_{2 n}=\sigma, \qquad & E_{1 t}-E_{2 t}=0 \\
B_{1 n}-B_{2 n}=0, \qquad & H_{1 t}-H_{2 t}=K.\end{array}
\end{equation}
Where we have fields $\b{E}_1, \b{D}_1, \b{B}_1, \b{H}_1$ in medium $1$ and fields $\b{E}_2, \b{D}_2, \b{B}_2, \b{H}_2$ in medium $2$. Subscript $t$ and $n$ denote the tangential and normal component respectively, $\sigma$ is the surface charge density and $\b K$ is the surface current density.

===Frequency dependence of the dielectric function===
As mentioned above the dielectric function $\varepsilon$ is a function of frequency. Some analytical properties of the dielectric function $\varepsilon(\omega)$ are well known, due to their relevance in optics, one example being the Kramers-Kronig relations. The dielectric function is in general, based on a nonlocal response function $\chi$, defined as
\begin{equation}
\varepsilon(\omega) = 1 + \int_{0}^{\infty} \chi(\tau) e^{i \omega \tau} \text{d}\tau.
\end{equation}
For a more detailed analysis of the analytical properties of the function, refer to any standard optics or advanced electromagnetics textbook. When it comes to numerical simulations, models of the dielectric function are of more interest than its analytical properties.
====Dielectric materials====
The simplest model for dielectrics is based on the equation of motion for the bound charges
\begin{equation}
m \frac{\mathrm{d}^{2} \mathbf{r}}{\mathrm{d} t^{2}}=-m \gamma \dot{\mathbf{r}}-m \omega_{0}^{2} \mathbf{r}+e \mathbf{E}(t)
\end{equation}
the terms on the right side describe the dissipation of energy, oscillations and the driving force of the external electric field. A short calculation leads to an expression of the frequency dependence of the dielectric function
\begin{equation}
\varepsilon(\omega)=1+ \frac{\omega_p}{\left(\omega_{0}^{2}-\omega^{2}\right)-i \gamma \omega}.
\end{equation}
The above expression can be simplified for high and low frequencies leading to the so-called Debye and Plasmonic relaxation models. The Lorentz model applies for frequencies where no terms in the denominator can be neglected. Typical Debye and Lorentz relaxations are shown in the plots below.

[[File:Debye.png|400px]][[File:Lorentz.png|400px]]

In general, the dielectric function of a material is modeled as a sum of many Debye and Lorentz relaxations, the parameters of which are determined by microscopic properties of the material itself, an extremely interesting case study is the dielectric function of water.
====Conductors====
When dealing with conductors we use the Drude model of conduction
\begin{equation}
m \frac{d \mathbf{v}(t)}{d t}=-m \gamma \mathbf{v}(t)+e \mathbf{E}(t)
\end{equation}
from which we quickly derive the equation for conductivity
\begin{equation}
\sigma(\omega)=\frac{n e^{2}}{m} \frac{1}{\gamma-i \omega},
\end{equation}
which we can relate to the dielectric function
\begin{equation}
\varepsilon(\omega)=1-\frac{\sigma(\omega)}{i \varepsilon_{0} \omega}.
\end{equation}
==Modeling in the time domain==
The single most frequent type of calculation used in computational electromagnetics is the simulation of electromagnetic fields in the time domain. The calculation propagates $\b E(\b r, t)$ and $\b H(\b r, t)$ in time, usually starting with all fields set to zero and some time dependent current source $\b J(\b r, t)$. The main advantage of time domain approaches is the ability to study transient phenomena while its main drawback is the study resonant problems.
===Finite differences in time domain (FDTD)===
==== Formulation ====
Finite differences in time domain were suggested in 1966 by Yee. The method uses finite differences to discretize Maxwell's equations and propagates the fields in time by alternating between updating the electric field strength and magnetic field strength. One of its defining attributes is the use of a staggered grid to represent $\b H$ and $\b E$ nodes. The staggered grid simplifies the calculation of curl expressions and helps with numerical stability and treating the boundaries of materials.

[[File:FDTD2.png|400px|center]]

This type of solution procedure provides an alternative to solving the wave equation for either $\b H$ or $\b E$.

==== Modeling sources ====
Since we are starting with all fields set to zero, we need to somehow introduce power into the model. This is usually done with current sources. The current signal used depends on the type of calculation. If we are interested in a steady-state solution we would use a sinusoidal source, that is slowly turned on and simulated the device until a steady state is reached. On the other hand, if we are interested in calculating the response of a device to many frequencies, we would use a Gaussian pulse and simulate the device until all the energy has left the system.

Taflove in his book [1] differentiates between five types of sources. The easiest to implement are so-called ''hard sources'', where one simply sets the value of a field in a node to the source value e.g.
\begin{equation}
\b E_z(i, t) = \sin 2\pi ft.
\end{equation}
The problem is that the node now acts as a perfect dielectric and incoming waves reflect from the point and do not pass through as they should. On the other hand one can use a ''soft source'' where energy is inserted into the domain via electric and magnetic (non-physical) current
\begin{equation}
\b E_z(i, t) = E_z(i, t) + J(t) \qquad \text{in} \qquad \b H_x(i, t) = H_x(i, t) + M(t).
\end{equation}
Waves no longer reflect from the source, however, the source still radiates in all directions, which is sometimes not desired. This issue is solved with TFSF (total field, scattered field) sources, which are based on the linearity of Maxwell's equations and allow to construct sources that inject the energy only in the desired direction. This is accomplished by decomposing the total fields into the scattered and incident
\begin{equation}
\b{E} = \b{E}^i + \b{E}^s \text{, } \qquad \b{H} = \b{H}^i + \b{H}^s.
\end{equation}
and then dividing the computational domain $\Omega$ into two subdomains. In one of the subdomains, we calculate only using the scattered part of the fields and in the other, we calculate with the total field. With careful treatment of the boundary between the subdomains one can construct an arbitrary incident wave. Alongside the already mentioned source types, ''waveguide'' sources and ''total scattered field'' sources. For an in-depth discussion refer to Taflove [1].

==== Transmittance and Reflectivity ====
Transmittance and reflectivity are dependent on the frequency of the incoming wave. One of the main strengths of time domain calculations is that one can calculate both coefficients for a wide frequency range with a single simulation. This is achieved by using a source that contains a wide range of frequencies i.e. a short Gaussian pulse.

Let us illustrate how to calculate $R$ and $T$ in one dimension. The calculation is similar in higher dimensions. We make two assumptions, firstly the source is directed directly into the device. This means that we will have only the reflected field on one side of the device and only the transmitted field on the other side of the device. Secondly, we assume that nothing reflects from the edges of the domain. We record the signals on both edges of the domain and label them $E_r(t)$ for reflected and $E_t(t)$ for transmitted. It is important to simulate for a long enough time that there is no energy left in the device. Otherwise, the calculation of the coefficients has not converged yet. The Fourier transforms of the signals at the edges, along with the transform of the source, define $R$ and $T$
\begin{equation}
R(\omega) = \left| \frac{E_o(\omega)}{S(\omega)} \right| \qquad \text{and} \qquad T(\omega)= \left| \frac{E_p(\omega)}{S(\omega)} \right|.
\end{equation}
Since we only know the field values at discrete time intervals the transformation is discrete
\begin{equation}
F(\omega) \cong \Delta t \sum_{m=1}^{M}\left(e^{-j \omega \Delta t}\right)^{m} \cdot f(m \Delta t),
\end{equation}
in practice, so that it is not needed to store the signals on both edges which may take up considerable computer memory, the sum can be calculated after each step in the simulation to save space.

<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}
</math>

==Modeling in the frequency domain==
In the frequency domain, two types of calculations are most common. Within the field of photonic crystals, bandstructures of periodic materials are sought by solving a generalized eigenvalue problem $Ax = \omega^2Bx$. In engineering applications, a response of a structure to a periodical source $\b J(x) e^{-i\omega t}$ is more often calculated. We will focus on the second type of calculation.

==== Formulation ====
We need to solve the equation that can be derived by combining Maxwell's equations and constitutive relations
\begin{equation}
\left[(\nabla \times \nabla \times)-\frac{\omega^{2}}{c^{2}} \varepsilon(\mathbf{r})\right] \mathbf{E}(\mathbf{r})=i \omega \mu_{0} \mathbf{J}(\mathbf{r}).
\end{equation}
The electric field strength $\b E(\b r)$ is complex, the time dependent solution is obtained via $\b E(\b r,t) = \b E(\b r) e^{i\omega t}$. The magnetic field can be obtained by taking a curl of the solution above
\begin{equation}
\mathbf{H}=-\frac{i}{\omega \mu_{0}} \nabla \times \mathbf{E}.
\end{equation}

==== Absorbing boundary conditions ====
Absorbing boundary conditions in the frequency domain can be realized in a number of ways. One of the most straight forward approaches is to impose the Sommerfeld radiation condition or a numerical approximation of it. However, it is more common to surround the computational domain with a thin absorbing layer, which does not reflect outgoing waves. The perfectly matched layer (PML) is one of such absorbing boundary conditions. There are many variants of PML's, but it is clear that the use of one worsens the conditioning of the problem matrix [2]. As it turns out UPML (Uniaxial PML) is numerically less stable than SC-PML (stretched coordinate PML), but efficient preconditioners exist for cases when SC-PML is hard to implement.

===== Uniaxial Perfectly Matched Layer =====
Inside the computational domain we have an amplitude equation for $\b E$,
\begin{equation}
\nabla \times \doubleunderline{\mu}^{-1} \nabla \times \mathbf{E}-\omega^{2} \varepsilon \mathbf{E}=-\mathrm{i} \omega \mathbf{J}.
\end{equation}
We surround the computational domain $\Omega$ with a thin domain $\Omega'$, where we have an uniaxial anisotropic material, which damps the electromagnetic field
\begin{equation}
\nabla \times \doubleunderline{\mu_{s}}^{-1} \nabla \times \mathbf{E}-\omega^{2} \doubleunderline{\varepsilon_{s}} \mathbf{E}=-\mathrm{i} \omega \mathbf{J}
\end{equation}
dielectric and permeability tensors are diagonal,
\begin{equation}
\doubleunderline{\varepsilon_{s}}=\varepsilon\left[\begin{array}{ccc}\frac{s_{y} s_{z}}{s_{x}} & 0 & 0 \\ 0 & \frac{s_{s} s_{x}}{s_{y}} & 0 \\ 0 & 0 & \frac{s_{x} s_{y}}{s_{z}}\end{array}\right], \quad \doubleunderline{\mu_{s}}=\mu\left[\begin{array}{ccc}\frac{s_{y} s_{z}}{s_{x}} & 0 & 0 \\ 0 & \frac{s_{z} s_{x}}{s_{y}} & 0 \\ 0 & 0 & \frac{s_{x} s_{y}}{s_{z}}\end{array}\right].
\end{equation}
To avoid reflections of waves from the boundary between layers, we slowly turn on the loss parameter $s_w$
\begin{equation}
s_{w}(l)=1-\mathrm{i} \frac{\sigma_{w}(l)}{\omega \varepsilon_{0}}
\end{equation}
where $l$ is the distance from the boundary of the PML.

===== Streched Coordinate Perfectly Matched Layer =====
In the stretched coordinate PML the coordinates are stretched and the PML no longer resembles any physical material. The equation inside the PML remains unchanged,
\begin{equation}
\label{eq:SC-PML}
\nabla_{s} \times \mu^{-1} \nabla_{s} \times \mathbf{E}-\omega^{2} \varepsilon \mathbf{E}=-\mathrm{i} \omega \mathbf{J}.
\end{equation}
However, the gradient is redefined as
\begin{equation}
\nabla_{s}=\hat{e}_x \frac{1}{s_{x}} \frac{\partial}{\partial x}+\hat{e}_y \frac{1}{s_{y}} \frac{\partial}{\partial y}+\hat{e}_z \frac{1}{s_{z}} \frac{\partial}{\partial z}.
\end{equation}
Just like with UPML, the damping is slowly turned on. SC-PML is considered more stable than UMPL, but it is usually trickier to implement. For a magnetically homogenous and isotropic medium the SC-PML equation simplifies to
\begin{equation}
-\nabla^2_s\b E + \nabla_s\left( \nabla_{s} \cdot \b E \right) -\omega^2 \eps \mu \b E = -i\omega\b J,
\end{equation}
in matrix form
\begin{multline}
-
\begin{bmatrix}
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_x}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_x}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_x}{z}\right) \\
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_y}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_y}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_y}{z}\right) \\
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_z}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_z}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_z}{z}\right)
\end{bmatrix}
+
\begin{bmatrix}
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_x}{x}\right) +
\frac{1}{s_y s_x}\frac{\partial E_y}{\partial x \partial y} +
\frac{1}{s_z s_x}\frac{\partial E_z}{\partial z \partial z} \\
\frac{1}{s_x s_y}\frac{\partial E_x}{\partial y \partial x} +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_y}{y}\right) +
\frac{1}{s_z s_y}\frac{\partial E_z}{\partial y \partial z} \\
\frac{1}{s_x s_z}\frac{\partial E_x}{\partial z \partial x} +
\frac{1}{s_y s_z}\frac{\partial E_y}{\partial z \partial y} +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_z}{z}\right)
\end{bmatrix} -
\\ -\omega^2 \eps \mu \b E
\begin{bmatrix}
E_x \\
E_y \\
E_z
\end{bmatrix}
= -i\omega
\begin{bmatrix}
J_x \\
J_y \\
J_z
\end{bmatrix}
\end{multline}
The expression further simplifies in cases where some of $s_w$'s are constant (cartesian coordinates).

[[File:radial_pml_2D_test_im.png|400px]][[File:radial_pml_2D_test_re.png|400px]][[File:radial_pml_2D_test_sfun.png|400px]]

The above graphs display a simple use case of SC-PML with radial geometry. Instead of being surrounded by a box, the computational domain is surrounded by an annulus of perfectly matched layer, that contains both $s_x$ and $s_y$ components and therefore dampens the wave in both $x$ and $y$ directions throughout the PML. Such approaches are not usually taken in conventional FDTD codes because of difficult implementation, this is not the case when using the Medusa library.

===== Benchmarking =====
Parameters of the PML are usually determined empirically, separately for each simulation, however, there are some rules of thumb. PMLs by Berenger's definition are exact for a continuous equation. When the equation is discretized reflections appear at the boundary of the PML. A lot of care needs to be taken so that the reflections are minimized. This is usually achieved by giving $\sigma_{w}(l)$ a polynomial relationship with the distance $l$ from the PML boundary
\begin{equation}
\sigma_{w}(l,) = \sigma_{max}\left(\frac{l}{d}\right)^m
\end{equation}

where $\sigma_{max}$ is the maximum loss parameter, $d$ is the thickness of the layer and $m$ is the power of the polynomial, and usually lies somewhere between $2$ and $4$. For a normally incident wave on the PML, we can choose the Reflectance of the PML and accordingly tune parameters $d$ and $m$ via
\begin{equation}
\sigma_{w, \max }=-\frac{(m+1) \ln R}{2 \eta_{0} d}
\end{equation}
where $\eta_0$ is the Impedance of free space.

==References==
# Taflove, A. and Hagness, S.C., 2005. Computational electrodynamics: the finite-difference time-domain method. Artech house.
# Shin, W. and Fan, S., 2012. Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell’s equations solvers. Journal of Computational Physics, 231(8), pp.3406-3431.
# Podgornik, Rudolf, and Andrej Vilfan. Elektromagnetno polje. DMFA-založništvo, 2012.
# John, S. G. J., et al. "Photonic crystals: molding the flow of light." Princeton University of press. 2008.

File:Radial pml 2D test im.png

2020-09-27T08:01:12Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:Radial pml 2D test re.png

2020-09-27T08:01:12Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:Radial pml 2D test sfun.png

2020-09-27T08:01:12Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

Point source near an anisotropic lens

2020-09-24T08:46:03Z

BlazS: /* Case */

Point source near an anisotropic lens

2020-09-24T08:43:33Z

BlazS: /* Case */

Scattering from an infinite cylinder

2020-09-24T08:19:59Z

BlazS: /* Case */

Back to [[Computational electromagnetics]].

==Case==
<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}</math>
As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. This case is important for two reasons. Firstly we will see how we can describe the solution in terms of scattered fields only, by including the incident wave in the boundary conditions. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
\begin{align}
\nabla^2 E_z^{int} +\eps_r \frac{\omega^2}{c_0^2}\thinspace E_z^{int} = 0 \qquad &\text{in} \quad D \label{eq:inner} \\
\nabla^2 E_z^s + \frac{\omega^2}{c_0^2}E_z^s = 0 \qquad & \text{outside} \quad D
\label{eq:outer}
\end{align}
The field $E_z^{int}$ is the total field inside, $E_z^s$ is only the scattered part of $E_z$ Outside the cylinder. The boundary conditions are
\begin{equation}
E_z^{int} - E_z^{s} =E_z^{inc} \qquad \text{on} \quad \partial D \label{eq:BC1}, \qquad \text{and} \qquad
\dpar{E_z^{int}}{n} - \frac{1}{\eps_r}\dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n} \qquad \text{on} \quad \partial D.
\end{equation}
The incident field $E_z^{inc}$ has an analytical form in our case a plane wave $e^{ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.

<syntaxhighlight lang="c++" line>
// Construct system
double l_in = epsr*std::pow(omega/c0, 2.0);
double l_out = std::pow(omega/c0, 2.0);
for (int i : cyl_dom.interior()){
op_inner.lap(i) + l_in*op_inner.value(i) = 0.0;
}
for (int i : outer_inter){
op_outer.lap(i) + l_out*op_outer.value(i) = 0.0;
}
for (int c = 0; c < interface_c_idx.size(); ++c){
int i = interface_c_idx[c]; // cyl index
int j = interface_o_idx[c]; // outer index
Vec2d pos = cyl_dom.pos(i);
double x = pos[0];
double y = pos[1];

// calculate normal derivative of the source funciton at the i-th node
// get normal, must point out of the cylinder
Vec2d normal = cyl_dom.normal(i);
std::complex<double> incident = std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
std::complex<double> din_dx =
1.0_i*k0*std::cos(t0)*std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
std::complex<double> din_dy =
1.0_i*k0*std::sin(t0)*std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
std::complex<double> dui_dn = normal[0]*din_dx + normal[1]*din_dy;

// continuity of the fields, and the incident field
op_inner.value(i) + (-1)*op_outer.value(j, Nouter + i) = incident;

// continuity of derivatives
op_outer.neumann(j, outer_dom.normal(j))
+ (1/epsr)*op_inner.neumann(i, cyl_dom.normal(i), j - Nouter)
= dui_dn; // todo, change this from Neumanns to derivatives!
}
// SC-PML - where either x or y is constant
for (int i : PML){
1.0/(sw[i]*sw[i])*op_outer.lap(i) + l_out*op_outer.value(i) +
((-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 0, i)*op_outer.der1(i, 0)
+(-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 1, i)*op_outer.der1(i, 1)) = 0.0;
}
for (int i: outer_bnd){
op_outer.value(i) = 0.0;
}

</syntaxhighlight>

==Analytical solution==
The solution to this problem is sought in terms of sums of cylindrical harmonics. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
\begin{equation}
\mathbf{E}_{i}=\widehat{z} E_{0} e^{-i k_{0} x}= \widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} J_{n}\left(k_{0} \rho\right) e^{i n \phi}.
\end{equation}
The scattered field is a sum of Hankel functions so it satisfies the Sommerfeld radiation condition
\begin{equation}
\mathbf{E}^{s}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} a_{n} H_{n}^{(2)}\left(k_{0} \rho\right) e^{i n \phi}, \quad \rho \geq a
\end{equation}
and finally the total field is again a sum ob Bessel functions of the first kind
\begin{equation}
\mathbf{E}^{t}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} b_{n} J_{n}\left(k_{1} \rho\right) e^{i n \phi}, \quad \rho \leq a.
\end{equation}
The boundary conditions for the field and its derivative on the boundary between the cylinder and free space surrounding it give two expressions
\begin{equation}
J_{n}\left(k_{0} a\right)+a_{n} H_{n}^{(2)}\left(k_{0} a\right)=b_{n} J_{n}\left(k_{1} a\right)
\end{equation}
and
\begin{equation}
k_0 J^{\prime}_{n}\left(k_{0} a\right)+ k_0 a_{n} H_{n}^{\prime(2)}\left(k_{0} a\right)=k_1 b_{n} J^{\prime}_{n}\left(k_{1} a\right).
\end{equation}
Which give the coefficients $a_n$ and $b_n$ in the above expansions as
\begin{equation}
a_{n}=-\frac{J_{n}\left(k_{0} a\right)-\Gamma_{h} J_{n}^{\prime}\left(k_{0} a\right)}{H_{n}^{(2)}\left(k_{0} a\right)-\Gamma_{h} H_{n}^{(2)^{\prime}}\left(k_{0} a\right)},
\end{equation}
where
\begin{equation}
\Gamma_{h}=\frac{Z_{0}}{Z_{1}} \frac{J_{n}\left(k_{1} a\right)}{J_{n}^{\prime}\left(k_{1} a\right)}
\end{equation}
and
\begin{equation}
b_n = \frac{J_n(k_0a)}{J_n(k_1a)} + a_n \frac{H_n^{(2)}(k_0a)}{J_n(k_1a)}.
\end{equation}

==Results==
[[File:FDFD_2D_dielcyl_ReEn_1.png|400px]][[File:FDFD_2D_dielcyl_ImEn_1.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_2.png|400px]][[File:FDFD_2D_dielcyl_ImEn_2.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_3.png|400px]][[File:FDFD_2D_dielcyl_ImEn_3.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_4.png|400px]][[File:FDFD_2D_dielcyl_ImEn_4.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_5.png|400px]][[File:FDFD_2D_dielcyl_ImEn_5.png|400px]]

[[File:FDFD_2D_benchmark_times_basic.png|800px]]

[[File:FDFD_2D_dx_conv2.png|800px]]

Scattering from an infinite cylinder

2020-09-24T08:17:52Z

BlazS: /* Case */

Back to [[Computational electromagnetics]].

==Case==
<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}</math>
As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. This case is important for two reasons. Firstly we will see how we can describe the solution in terms of scattered fields only, by including the incident wave in the boundary conditions. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
\begin{align}
\nabla^2 E_z^{int} +\eps_r \frac{\omega^2}{c_0^2}\thinspace E_z^{int} = 0 \qquad &\text{in} \quad D \label{eq:inner} \\
\nabla^2 E_z^s + \frac{\omega^2}{c_0^2}E_z^s = 0 \qquad & \text{outside} \quad D
\label{eq:outer}
\end{align}
The field $E_z^{int}$ is the total field inside, $E_z^s$ is only the scattered part of $E_z$ Outside the cylinder. The boundary conditions are
\begin{equation}
E_z^{int} - E_z^{s} =E_z^{inc} \qquad \text{on} \quad \partial D \label{eq:BC1}, \qquad \text{and} \qquad
\dpar{E_z^{int}}{n} - \frac{1}{\eps_r}\dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n} \qquad \text{on} \quad \partial D.
\end{equation}
The incident field $E_z^{inc}$ has an analytical form in our case a plane wave $e^{ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.

<syntaxhighlight lang="c++" line>
// Construct system
double l_in = epsr*std::pow(omega/c0, 2.0);
double l_out = std::pow(omega/c0, 2.0);
for (int i : cyl_dom.interior()){
op_inner.lap(i) + l_in*op_inner.value(i) = 0.0;
}
for (int i : outer_inter){
op_outer.lap(i) + l_out*op_outer.value(i) = 0.0;
}
for (int c = 0; c < interface_c_idx.size(); ++c){
int i = interface_c_idx[c]; // cyl index
int j = interface_o_idx[c]; // outer index
Vec2d pos = cyl_dom.pos(i);
double x = pos[0];
double y = pos[1];

// calculate normal derivative of the source funciton at the i-th node
// get normal, must point out of the cylinder
Vec2d normal = cyl_dom.normal(i);
std::complex<double> incident = std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
std::complex<double> din_dx = 1.0_i*k0*std::cos(t0)*std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
std::complex<double> din_dy = 1.0_i*k0*std::sin(t0)*std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
std::complex<double> dui_dn = normal[0]*din_dx + normal[1]*din_dy;

// continuity of the fields, and the incident field
op_inner.value(i) + (-1)*op_outer.value(j, Nouter + i) = incident;

// continuity of derivatives
op_outer.neumann(j, outer_dom.normal(j)) + (1/epsr)*op_inner.neumann(i, cyl_dom.normal(i), j - Nouter)
= dui_dn;
}
// SC-PML - where either x or y is constant
for (int i : PML){
1.0/(sw[i]*sw[i])*op_outer.lap(i) + l_out*op_outer.value(i) +
((-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 0, i)*op_outer.der1(i, 0)
+(-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 1, i)*op_outer.der1(i, 1)) = 0.0;
}
for (int i: outer_bnd){
op_outer.value(i) = 0.0;
}

</syntaxhighlight>

==Analytical solution==
The solution to this problem is sought in terms of sums of cylindrical harmonics. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
\begin{equation}
\mathbf{E}_{i}=\widehat{z} E_{0} e^{-i k_{0} x}= \widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} J_{n}\left(k_{0} \rho\right) e^{i n \phi}.
\end{equation}
The scattered field is a sum of Hankel functions so it satisfies the Sommerfeld radiation condition
\begin{equation}
\mathbf{E}^{s}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} a_{n} H_{n}^{(2)}\left(k_{0} \rho\right) e^{i n \phi}, \quad \rho \geq a
\end{equation}
and finally the total field is again a sum ob Bessel functions of the first kind
\begin{equation}
\mathbf{E}^{t}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} b_{n} J_{n}\left(k_{1} \rho\right) e^{i n \phi}, \quad \rho \leq a.
\end{equation}
The boundary conditions for the field and its derivative on the boundary between the cylinder and free space surrounding it give two expressions
\begin{equation}
J_{n}\left(k_{0} a\right)+a_{n} H_{n}^{(2)}\left(k_{0} a\right)=b_{n} J_{n}\left(k_{1} a\right)
\end{equation}
and
\begin{equation}
k_0 J^{\prime}_{n}\left(k_{0} a\right)+ k_0 a_{n} H_{n}^{\prime(2)}\left(k_{0} a\right)=k_1 b_{n} J^{\prime}_{n}\left(k_{1} a\right).
\end{equation}
Which give the coefficients $a_n$ and $b_n$ in the above expansions as
\begin{equation}
a_{n}=-\frac{J_{n}\left(k_{0} a\right)-\Gamma_{h} J_{n}^{\prime}\left(k_{0} a\right)}{H_{n}^{(2)}\left(k_{0} a\right)-\Gamma_{h} H_{n}^{(2)^{\prime}}\left(k_{0} a\right)},
\end{equation}
where
\begin{equation}
\Gamma_{h}=\frac{Z_{0}}{Z_{1}} \frac{J_{n}\left(k_{1} a\right)}{J_{n}^{\prime}\left(k_{1} a\right)}
\end{equation}
and
\begin{equation}
b_n = \frac{J_n(k_0a)}{J_n(k_1a)} + a_n \frac{H_n^{(2)}(k_0a)}{J_n(k_1a)}.
\end{equation}

==Results==
[[File:FDFD_2D_dielcyl_ReEn_1.png|400px]][[File:FDFD_2D_dielcyl_ImEn_1.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_2.png|400px]][[File:FDFD_2D_dielcyl_ImEn_2.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_3.png|400px]][[File:FDFD_2D_dielcyl_ImEn_3.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_4.png|400px]][[File:FDFD_2D_dielcyl_ImEn_4.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_5.png|400px]][[File:FDFD_2D_dielcyl_ImEn_5.png|400px]]

[[File:FDFD_2D_benchmark_times_basic.png|800px]]

[[File:FDFD_2D_dx_conv2.png|800px]]

Scattering from an infinite cylinder

2020-09-24T07:54:00Z

BlazS: /* Case */

Back to [[Computational electromagnetics]].

==Case==
<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}</math>
As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. This case is important for two reasons. Firstly we will see how we can describe the solution in terms of scattered fields only, by including the incident wave in the boundary conditions. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
\begin{align}
\nabla^2 E_z^{int} +\eps_r \frac{\omega^2}{c_0^2}\thinspace E_z^{int} = 0 \qquad &\text{in} \quad D \label{eq:inner} \\
\nabla^2 E_z^s + \frac{\omega^2}{c_0^2}E_z^s = 0 \qquad & \text{outside} \quad D
\label{eq:outer}
\end{align}
The field $E_z^{int}$ is the total field inside, $E_z^s$ is only the scattered part of $E_z$ Outside the cylinder. The boundary conditions are
\begin{equation}
E_z^{int} - E_z^{s} =E_z^{inc} \qquad \text{on} \quad \partial D \label{eq:BC1}, \qquad \text{and} \qquad
\dpar{E_z^{int}}{n} - \frac{1}{\eps_r}\dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n} \qquad \text{on} \quad \partial D.
\end{equation}
The incident field $E_z^{inc}$ has an analytical form in our case a plane wave $e^{ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.

==Analytical solution==
The solution to this problem is sought in terms of sums of cylindrical harmonics. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
\begin{equation}
\mathbf{E}_{i}=\widehat{z} E_{0} e^{-i k_{0} x}= \widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} J_{n}\left(k_{0} \rho\right) e^{i n \phi}.
\end{equation}
The scattered field is a sum of Hankel functions so it satisfies the Sommerfeld radiation condition
\begin{equation}
\mathbf{E}^{s}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} a_{n} H_{n}^{(2)}\left(k_{0} \rho\right) e^{i n \phi}, \quad \rho \geq a
\end{equation}
and finally the total field is again a sum ob Bessel functions of the first kind
\begin{equation}
\mathbf{E}^{t}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} b_{n} J_{n}\left(k_{1} \rho\right) e^{i n \phi}, \quad \rho \leq a.
\end{equation}
The boundary conditions for the field and its derivative on the boundary between the cylinder and free space surrounding it give two expressions
\begin{equation}
J_{n}\left(k_{0} a\right)+a_{n} H_{n}^{(2)}\left(k_{0} a\right)=b_{n} J_{n}\left(k_{1} a\right)
\end{equation}
and
\begin{equation}
k_0 J^{\prime}_{n}\left(k_{0} a\right)+ k_0 a_{n} H_{n}^{\prime(2)}\left(k_{0} a\right)=k_1 b_{n} J^{\prime}_{n}\left(k_{1} a\right).
\end{equation}
Which give the coefficients $a_n$ and $b_n$ in the above expansions as
\begin{equation}
a_{n}=-\frac{J_{n}\left(k_{0} a\right)-\Gamma_{h} J_{n}^{\prime}\left(k_{0} a\right)}{H_{n}^{(2)}\left(k_{0} a\right)-\Gamma_{h} H_{n}^{(2)^{\prime}}\left(k_{0} a\right)},
\end{equation}
where
\begin{equation}
\Gamma_{h}=\frac{Z_{0}}{Z_{1}} \frac{J_{n}\left(k_{1} a\right)}{J_{n}^{\prime}\left(k_{1} a\right)}
\end{equation}
and
\begin{equation}
b_n = \frac{J_n(k_0a)}{J_n(k_1a)} + a_n \frac{H_n^{(2)}(k_0a)}{J_n(k_1a)}.
\end{equation}

==Results==
[[File:FDFD_2D_dielcyl_ReEn_1.png|400px]][[File:FDFD_2D_dielcyl_ImEn_1.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_2.png|400px]][[File:FDFD_2D_dielcyl_ImEn_2.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_3.png|400px]][[File:FDFD_2D_dielcyl_ImEn_3.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_4.png|400px]][[File:FDFD_2D_dielcyl_ImEn_4.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_5.png|400px]][[File:FDFD_2D_dielcyl_ImEn_5.png|400px]]

[[File:FDFD_2D_benchmark_times_basic.png|800px]]

[[File:FDFD_2D_dx_conv2.png|800px]]

Scattering from an infinite cylinder

2020-09-24T07:53:06Z

BlazS:

Back to [[Computational electromagnetics]].

==Case==
<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}</math>
As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. This case is important for two reasons. Firstly we will see how we can describe the solution in terms of scattered fields only, by including the incident wave in the boundary conditions. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
\begin{align}
\nabla^2 E_z^{int} +\eps_r \frac{\omega^2}{c_0^2}\thinspace E_z^{int} = 0 \qquad &\text{in} \quad D \label{eq:inner} \\
\nabla^2 E_z^s + \frac{\omega^2}{c_0^2}E_z^s = 0 \qquad & \text{outside} \quad D
\label{eq:outer}
\end{align}
The field $E_z^{int}$ is the total field inside, $E_z^s$ is only the scattered part of $E_z$ Outside the cylinder. The boundary conditions are
\begin{equation}
E_z^{int} - E_z^{s} =E_z^{inc} \qquad \text{on} \quad \partial \Omega \label{eq:BC1}, \qquad \text{and} \qquad
\dpar{E_z^{int}}{n} - \frac{1}{\eps_r}\dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n} \qquad \text{on} \quad \partial \Omega.
\end{equation}
The incident field $u^i$ has an analytical form in our case a plane wave $e^{ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.

==Analytical solution==
The solution to this problem is sought in terms of sums of cylindrical harmonics. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
\begin{equation}
\mathbf{E}_{i}=\widehat{z} E_{0} e^{-i k_{0} x}= \widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} J_{n}\left(k_{0} \rho\right) e^{i n \phi}.
\end{equation}
The scattered field is a sum of Hankel functions so it satisfies the Sommerfeld radiation condition
\begin{equation}
\mathbf{E}^{s}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} a_{n} H_{n}^{(2)}\left(k_{0} \rho\right) e^{i n \phi}, \quad \rho \geq a
\end{equation}
and finally the total field is again a sum ob Bessel functions of the first kind
\begin{equation}
\mathbf{E}^{t}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} b_{n} J_{n}\left(k_{1} \rho\right) e^{i n \phi}, \quad \rho \leq a.
\end{equation}
The boundary conditions for the field and its derivative on the boundary between the cylinder and free space surrounding it give two expressions
\begin{equation}
J_{n}\left(k_{0} a\right)+a_{n} H_{n}^{(2)}\left(k_{0} a\right)=b_{n} J_{n}\left(k_{1} a\right)
\end{equation}
and
\begin{equation}
k_0 J^{\prime}_{n}\left(k_{0} a\right)+ k_0 a_{n} H_{n}^{\prime(2)}\left(k_{0} a\right)=k_1 b_{n} J^{\prime}_{n}\left(k_{1} a\right).
\end{equation}
Which give the coefficients $a_n$ and $b_n$ in the above expansions as
\begin{equation}
a_{n}=-\frac{J_{n}\left(k_{0} a\right)-\Gamma_{h} J_{n}^{\prime}\left(k_{0} a\right)}{H_{n}^{(2)}\left(k_{0} a\right)-\Gamma_{h} H_{n}^{(2)^{\prime}}\left(k_{0} a\right)},
\end{equation}
where
\begin{equation}
\Gamma_{h}=\frac{Z_{0}}{Z_{1}} \frac{J_{n}\left(k_{1} a\right)}{J_{n}^{\prime}\left(k_{1} a\right)}
\end{equation}
and
\begin{equation}
b_n = \frac{J_n(k_0a)}{J_n(k_1a)} + a_n \frac{H_n^{(2)}(k_0a)}{J_n(k_1a)}.
\end{equation}

==Results==
[[File:FDFD_2D_dielcyl_ReEn_1.png|400px]][[File:FDFD_2D_dielcyl_ImEn_1.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_2.png|400px]][[File:FDFD_2D_dielcyl_ImEn_2.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_3.png|400px]][[File:FDFD_2D_dielcyl_ImEn_3.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_4.png|400px]][[File:FDFD_2D_dielcyl_ImEn_4.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_5.png|400px]][[File:FDFD_2D_dielcyl_ImEn_5.png|400px]]

[[File:FDFD_2D_benchmark_times_basic.png|800px]]

[[File:FDFD_2D_dx_conv2.png|800px]]

Scattering from an infinite cylinder

2020-09-24T07:40:25Z

BlazS:

Back to [[Computational electromagnetics]].

==Case==
<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}</math>
As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. This case is important for two reasons. Firstly we will see how we can describe the solution in terms of scattered fields only, by including the incident wave in the boundary conditions. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
\begin{align}
\nabla^2 E_z^{int} +\eps_r \frac{\omega^2}{c_0^2}\thinspace E_z^{int} = 0 \qquad &\text{and} \quad \Omega \label{eq:inner} \\
\nabla^2 E_z^s + \frac{\omega^2}{c_0^2}E_z^s = 0 \qquad & \text{outside} \quad \Omega
\label{eq:outer}
\end{align}
The field $v$ is the total field inside, $u^s$ is only the scattered part of $E_z$ Outside the cylinder. The boundary conditions are
\begin{equation}
E_z^{int} - E_z^{s} =u^{inc} \qquad \text{on} \quad \partial \Omega \label{eq:BC1}, \qquad \text{and} \qquad
\dpar{E_z^{int}}{n} - \frac{1}{\eps_r}\dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n} \qquad \text{on} \quad \partial \Omega.
\end{equation}
The incident field $u^i$ has an analytical form in our case a plane wave $e^{ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.

==Analytical solution==
The solution to this problem is sought in terms of sums of cylindrical harmonics. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
\begin{equation}
\mathbf{E}_{i}=\widehat{z} E_{0} e^{-i k_{0} x}= \widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} J_{n}\left(k_{0} \rho\right) e^{i n \phi}.
\end{equation}
The scattered field is a sum of Hankel functions so it satisfies the Sommerfeld radiation condition
\begin{equation}
\mathbf{E}^{s}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} a_{n} H_{n}^{(2)}\left(k_{0} \rho\right) e^{i n \phi}, \quad \rho \geq a
\end{equation}
and finally the total field is again a sum ob Bessel functions of the first kind
\begin{equation}
\mathbf{E}^{t}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} b_{n} J_{n}\left(k_{1} \rho\right) e^{i n \phi}, \quad \rho \leq a.
\end{equation}
The boundary conditions for the field and its derivative on the boundary between the cylinder and free space surrounding it give two expressions
\begin{equation}
J_{n}\left(k_{0} a\right)+a_{n} H_{n}^{(2)}\left(k_{0} a\right)=b_{n} J_{n}\left(k_{1} a\right)
\end{equation}
and
\begin{equation}
k_0 J^{\prime}_{n}\left(k_{0} a\right)+ k_0 a_{n} H_{n}^{\prime(2)}\left(k_{0} a\right)=k_1 b_{n} J^{\prime}_{n}\left(k_{1} a\right).
\end{equation}
Which give the coefficients $a_n$ and $b_n$ in the above expansions as
\begin{equation}
a_{n}=-\frac{J_{n}\left(k_{0} a\right)-\Gamma_{h} J_{n}^{\prime}\left(k_{0} a\right)}{H_{n}^{(2)}\left(k_{0} a\right)-\Gamma_{h} H_{n}^{(2)^{\prime}}\left(k_{0} a\right)},
\end{equation}
where
\begin{equation}
\Gamma_{h}=\frac{Z_{0}}{Z_{1}} \frac{J_{n}\left(k_{1} a\right)}{J_{n}^{\prime}\left(k_{1} a\right)}
\end{equation}
and
\begin{equation}
b_n = \frac{J_n(k_0a)}{J_n(k_1a)} + a_n \frac{H_n^{(2)}(k_0a)}{J_n(k_1a)}.
\end{equation}

==Results==
[[File:FDFD_2D_dielcyl_ReEn_1.png|400px]][[File:FDFD_2D_dielcyl_ImEn_1.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_2.png|400px]][[File:FDFD_2D_dielcyl_ImEn_2.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_3.png|400px]][[File:FDFD_2D_dielcyl_ImEn_3.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_4.png|400px]][[File:FDFD_2D_dielcyl_ImEn_4.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_5.png|400px]][[File:FDFD_2D_dielcyl_ImEn_5.png|400px]]

[[File:FDFD_2D_benchmark_times_basic.png|800px]]

[[File:FDFD_2D_dx_conv2.png|800px]]

Scattering from an infinite cylinder

2020-09-24T07:38:23Z

BlazS:

Back to [[Computational electromagnetics]].

==Case==
<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}</math>
As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. Zgled je pomemben iz dveh razlogov, prvi\v c, ogledali si bomo kako vpadno polje zapi\v semo v rob med materialom in prostim prostorom in lahko kot take vire analiti\v cno predstavimo\footnote{Tu to ne pride posebno prav, ker imamo opravka z ravnim valom, v primerih ko imamo to\v ckaste izvire pa je zelo prakti\v cno.}. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
\begin{align}
\nabla^2 E_z^{int} +\eps_r \frac{\omega^2}{c_0^2}\thinspace E_z^{int} = 0 \qquad &\text{and} \quad \Omega \label{eq:inner} \\
\nabla^2 E_z^s + \frac{\omega^2}{c_0^2}E_z^s = 0 \qquad & \text{outside} \quad \Omega
\label{eq:outer}
\end{align}
The field $v$ is the total field inside, $u^s$ is only the scattered part of $E_z$ Outside the cylinder. The boundary conditions are
\begin{equation}
E_z^{int} - E_z^{s} =u^{inc} \qquad \text{on} \quad \partial \Omega \label{eq:BC1}, \qquad \text{and} \qquad
\dpar{E_z^{int}}{n} - \frac{1}{\eps_r}\dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n} \qquad \text{on} \quad \partial \Omega.
\end{equation}
The incident field $u^i$ has an analytical form in our case a plane wave $e^{ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.

==Analytical solution==
The solution to this problem is sought in terms of sums of cylindrical harmonics. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
\begin{equation}
\mathbf{E}_{i}=\widehat{z} E_{0} e^{-i k_{0} x}= \widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} J_{n}\left(k_{0} \rho\right) e^{i n \phi}.
\end{equation}
The scattered field is a sum of Hankel functions so it satisfies the Sommerfeld radiation condition
\begin{equation}
\mathbf{E}^{s}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} a_{n} H_{n}^{(2)}\left(k_{0} \rho\right) e^{i n \phi}, \quad \rho \geq a
\end{equation}
and finally the total field is again a sum ob Bessel functions of the first kind
\begin{equation}
\mathbf{E}^{t}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} b_{n} J_{n}\left(k_{1} \rho\right) e^{i n \phi}, \quad \rho \leq a.
\end{equation}
The boundary conditions for the field and its derivative on the boundary between the cylinder and free space surrounding it give two expressions
\begin{equation}
J_{n}\left(k_{0} a\right)+a_{n} H_{n}^{(2)}\left(k_{0} a\right)=b_{n} J_{n}\left(k_{1} a\right)
\end{equation}
and
\begin{equation}
k_0 J^{\prime}_{n}\left(k_{0} a\right)+ k_0 a_{n} H_{n}^{\prime(2)}\left(k_{0} a\right)=k_1 b_{n} J^{\prime}_{n}\left(k_{1} a\right).
\end{equation}
Which give the coefficients $a_n$ and $b_n$ in the above expansions as
\begin{equation}
a_{n}=-\frac{J_{n}\left(k_{0} a\right)-\Gamma_{h} J_{n}^{\prime}\left(k_{0} a\right)}{H_{n}^{(2)}\left(k_{0} a\right)-\Gamma_{h} H_{n}^{(2)^{\prime}}\left(k_{0} a\right)},
\end{equation}
where
\begin{equation}
\Gamma_{h}=\frac{Z_{0}}{Z_{1}} \frac{J_{n}\left(k_{1} a\right)}{J_{n}^{\prime}\left(k_{1} a\right)}
\end{equation}
and
\begin{equation}
b_n = \frac{J_n(k_0a)}{J_n(k_1a)} + a_n \frac{H_n^{(2)}(k_0a)}{J_n(k_1a)}.
\end{equation}

==Results==
[[File:FDFD_2D_dielcyl_ReEn_1.png|400px]][[File:FDFD_2D_dielcyl_ImEn_1.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_2.png|400px]][[File:FDFD_2D_dielcyl_ImEn_2.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_3.png|400px]][[File:FDFD_2D_dielcyl_ImEn_3.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_4.png|400px]][[File:FDFD_2D_dielcyl_ImEn_4.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_5.png|400px]][[File:FDFD_2D_dielcyl_ImEn_5.png|400px]]

[[File:FDFD_2D_benchmark_times_basic.png|800px]]

[[File:FDFD_2D_dx_conv2.png|800px]]

Triple dielectric step in 1D

2020-08-14T09:04:29Z

BlazS: /* Results */

Back to [[Computational electromagnetics]].
== Case ==
<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
</math>
A good starting point in the study of computational electromagnetics is a calculation of reflectance and transmissivity of a simple device. In this example, we will be studying device made of three dielectric layers with different dielectric constants. The layers are infinite and thus the problem simplifies to one dimension. We still need to differentiate between two different cases, the incident wave can be either $y$ polarized
\begin{equation}
\dpar{H_z}{t} =- \frac{1}{\mu(r)}\dpar{E_y}{x}, \qquad \text{in} \qquad \dpar{E_y}{t} = -\frac{1}{\eps(r)}\dpar{H_z}{x},
\end{equation}
or $z$ polarized
\begin{equation}
\dpar{H_y}{t} = \frac{1}{\mu(r)}\dpar{E_z}{x}, \qquad \text{in} \qquad \dpar{E_y}{t} = \frac{1}{\eps(r)}\dpar{H_y}{x},
\end{equation}

== Results ==
The problem for $z$-polarized incident light was solved with the medusa library using $N = 4001$ points and $\Delta t = dx/2c_0$. The simulation used a TFSF source that was located at $x = 0.1$, a Gaussian source that contained frequencies up to 2.5GHz was used. The plots below show a single step during the simulation and the transmissivity and reflectance of the device.

[[File:eps_step1d-mainEn.png|800px]]

[[File:eps_step1dEn-TR.png|800px]]

File:Eps step1dEn-TR.png

2020-08-14T09:04:19Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:Eps step1d-mainEn.png

2020-08-14T09:04:18Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

Scattering from an infinite cylinder

2020-08-14T07:26:47Z

BlazS: /* Results */

Back to [[Computational electromagnetics]].

==Case==
<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}</math>
As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. Zgled je pomemben iz dveh razlogov, prvi\v c, ogledali si bomo kako vpadno polje zapi\v semo v rob med materialom in prostim prostorom in lahko kot take vire analiti\v cno predstavimo\footnote{Tu to ne pride posebno prav, ker imamo opravka z ravnim valom, v primerih ko imamo to\v ckaste izvire pa je zelo prakti\v cno.}. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
\begin{align}
\nabla^2 v +\eps_r \frac{\omega^2}{c_0^2}\thinspace v = 0 \qquad &\text{and} \quad \Omega \label{eq:inner} \\
\nabla^2 u^s + \frac{\omega^2}{c_0^2}u^s = 0 \qquad & \text{outside} \quad \Omega
\label{eq:outer}
\end{align}
The field $v$ is the total field inside, $u^s$ is only the scattered part of $E_z$ Outside the cylinder. The boundary conditions are
\begin{equation}
v - u^s =u^i \qquad \text{on} \quad \partial \Omega \label{eq:BC1}, \qquad \text{and} \qquad
\dpar{v}{n} - \frac{1}{\eps_r}\dpar{u^s}{n} = \dpar{u^i}{n} \qquad \text{on} \quad \partial \Omega.
\end{equation}
The incident field $u^i$ has an analytical form in our case a plane wave $e^{ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.

==Analytical solution==
The solution to this problem is sought in terms of sums of cylindrical harmonics. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
\begin{equation}
\mathbf{E}_{i}=\widehat{z} E_{0} e^{-i k_{0} x}= \widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} J_{n}\left(k_{0} \rho\right) e^{i n \phi}.
\end{equation}
The scattered field is a sum of Hankel functions so it satisfies the Sommerfeld radiation condition
\begin{equation}
\mathbf{E}^{s}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} a_{n} H_{n}^{(2)}\left(k_{0} \rho\right) e^{i n \phi}, \quad \rho \geq a
\end{equation}
and finally the total field is again a sum ob Bessel functions of the first kind
\begin{equation}
\mathbf{E}^{t}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} b_{n} J_{n}\left(k_{1} \rho\right) e^{i n \phi}, \quad \rho \leq a.
\end{equation}
The boundary conditions for the field and its derivative on the boundary between the cylinder and free space surrounding it give two expressions
\begin{equation}
J_{n}\left(k_{0} a\right)+a_{n} H_{n}^{(2)}\left(k_{0} a\right)=b_{n} J_{n}\left(k_{1} a\right)
\end{equation}
and
\begin{equation}
k_0 J^{\prime}_{n}\left(k_{0} a\right)+ k_0 a_{n} H_{n}^{\prime(2)}\left(k_{0} a\right)=k_1 b_{n} J^{\prime}_{n}\left(k_{1} a\right).
\end{equation}
Which give the coefficients $a_n$ and $b_n$ in the above expansions as
\begin{equation}
a_{n}=-\frac{J_{n}\left(k_{0} a\right)-\Gamma_{h} J_{n}^{\prime}\left(k_{0} a\right)}{H_{n}^{(2)}\left(k_{0} a\right)-\Gamma_{h} H_{n}^{(2)^{\prime}}\left(k_{0} a\right)},
\end{equation}
where
\begin{equation}
\Gamma_{h}=\frac{Z_{0}}{Z_{1}} \frac{J_{n}\left(k_{1} a\right)}{J_{n}^{\prime}\left(k_{1} a\right)}
\end{equation}
and
\begin{equation}
b_n = \frac{J_n(k_0a)}{J_n(k_1a)} + a_n \frac{H_n^{(2)}(k_0a)}{J_n(k_1a)}.
\end{equation}

==Results==
[[File:FDFD_2D_dielcyl_ReEn_1.png|400px]][[File:FDFD_2D_dielcyl_ImEn_1.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_2.png|400px]][[File:FDFD_2D_dielcyl_ImEn_2.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_3.png|400px]][[File:FDFD_2D_dielcyl_ImEn_3.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_4.png|400px]][[File:FDFD_2D_dielcyl_ImEn_4.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_5.png|400px]][[File:FDFD_2D_dielcyl_ImEn_5.png|400px]]

[[File:FDFD_2D_benchmark_times_basic.png|800px]]

[[File:FDFD_2D_dx_conv2.png|800px]]

File:FDFD 2D dielcyl ReEn 5.png

2020-08-14T07:25:07Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:FDFD 2D dielcyl ReEn 3.png

2020-08-14T07:25:06Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:FDFD 2D dielcyl ReEn 4.png

2020-08-14T07:25:06Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:FDFD 2D dielcyl ReEn 1.png

2020-08-14T07:25:05Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:FDFD 2D dielcyl ReEn 2.png

2020-08-14T07:25:05Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:FDFD 2D dielcyl ImEn 4.png

2020-08-14T07:25:04Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:FDFD 2D dielcyl ImEn 5.png

2020-08-14T07:25:04Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:FDFD 2D dielcyl ImEn 1.png

2020-08-14T07:25:03Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:FDFD 2D dielcyl ImEn 2.png

2020-08-14T07:25:03Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:FDFD 2D dielcyl ImEn 3.png

2020-08-14T07:25:03Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

Point source near an anisotropic lens

2020-08-14T07:06:22Z

BlazS: /* Results */

File:FDFD 2D anizolens ReEn.png

2020-08-14T07:06:13Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:FDFD 2D anizolens ImEn.png

2020-08-14T07:05:44Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:FDFD 2D anizolens Im.png

2020-08-14T07:04:58Z

BlazS: BlazS uploaded a new version of File:FDFD 2D anizolens Im.png

File uploaded with MsUpload

File:FDFD 2D anizolens Re.png

2020-08-14T07:03:48Z

BlazS: BlazS uploaded a new version of File:FDFD 2D anizolens Re.png

File uploaded with MsUpload

File:FDFD 2D anizolens Im.png

2020-08-14T07:03:47Z

BlazS: BlazS uploaded a new version of File:FDFD 2D anizolens Im.png

File uploaded with MsUpload

Scattering from an infinite cylinder

2020-08-11T17:56:35Z

BlazS: /* Analytical solution */

Back to [[Computational electromagnetics]].

==Case==
<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}</math>
As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. Zgled je pomemben iz dveh razlogov, prvi\v c, ogledali si bomo kako vpadno polje zapi\v semo v rob med materialom in prostim prostorom in lahko kot take vire analiti\v cno predstavimo\footnote{Tu to ne pride posebno prav, ker imamo opravka z ravnim valom, v primerih ko imamo to\v ckaste izvire pa je zelo prakti\v cno.}. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
\begin{align}
\nabla^2 v +\eps_r \frac{\omega^2}{c_0^2}\thinspace v = 0 \qquad &\text{and} \quad \Omega \label{eq:inner} \\
\nabla^2 u^s + \frac{\omega^2}{c_0^2}u^s = 0 \qquad & \text{outside} \quad \Omega
\label{eq:outer}
\end{align}
The field $v$ is the total field inside, $u^s$ is only the scattered part of $E_z$ Outside the cylinder. The boundary conditions are
\begin{equation}
v - u^s =u^i \qquad \text{on} \quad \partial \Omega \label{eq:BC1}, \qquad \text{and} \qquad
\dpar{v}{n} - \frac{1}{\eps_r}\dpar{u^s}{n} = \dpar{u^i}{n} \qquad \text{on} \quad \partial \Omega.
\end{equation}
The incident field $u^i$ has an analytical form in our case a plane wave $e^{ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.

==Analytical solution==
The solution to this problem is sought in terms of sums of cylindrical harmonics. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
\begin{equation}
\mathbf{E}_{i}=\widehat{z} E_{0} e^{-i k_{0} x}= \widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} J_{n}\left(k_{0} \rho\right) e^{i n \phi}.
\end{equation}
The scattered field is a sum of Hankel functions so it satisfies the Sommerfeld radiation condition
\begin{equation}
\mathbf{E}^{s}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} a_{n} H_{n}^{(2)}\left(k_{0} \rho\right) e^{i n \phi}, \quad \rho \geq a
\end{equation}
and finally the total field is again a sum ob Bessel functions of the first kind
\begin{equation}
\mathbf{E}^{t}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} b_{n} J_{n}\left(k_{1} \rho\right) e^{i n \phi}, \quad \rho \leq a.
\end{equation}
The boundary conditions for the field and its derivative on the boundary between the cylinder and free space surrounding it give two expressions
\begin{equation}
J_{n}\left(k_{0} a\right)+a_{n} H_{n}^{(2)}\left(k_{0} a\right)=b_{n} J_{n}\left(k_{1} a\right)
\end{equation}
and
\begin{equation}
k_0 J^{\prime}_{n}\left(k_{0} a\right)+ k_0 a_{n} H_{n}^{\prime(2)}\left(k_{0} a\right)=k_1 b_{n} J^{\prime}_{n}\left(k_{1} a\right).
\end{equation}
Which give the coefficients $a_n$ and $b_n$ in the above expansions as
\begin{equation}
a_{n}=-\frac{J_{n}\left(k_{0} a\right)-\Gamma_{h} J_{n}^{\prime}\left(k_{0} a\right)}{H_{n}^{(2)}\left(k_{0} a\right)-\Gamma_{h} H_{n}^{(2)^{\prime}}\left(k_{0} a\right)},
\end{equation}
where
\begin{equation}
\Gamma_{h}=\frac{Z_{0}}{Z_{1}} \frac{J_{n}\left(k_{1} a\right)}{J_{n}^{\prime}\left(k_{1} a\right)}
\end{equation}
and
\begin{equation}
b_n = \frac{J_n(k_0a)}{J_n(k_1a)} + a_n \frac{H_n^{(2)}(k_0a)}{J_n(k_1a)}.
\end{equation}

==Results==
[[File:FDFD_2D_dielcyl_Re_1.png|400px]][[File:FDFD_2D_dielcyl_Im_1.png|400px]]

[[File:FDFD_2D_dielcyl_Re_2.png|400px]][[File:FDFD_2D_dielcyl_Im_2.png|400px]]

[[File:FDFD_2D_dielcyl_Re_3.png|400px]][[File:FDFD_2D_dielcyl_Im_3.png|400px]]

[[File:FDFD_2D_dielcyl_Re_4.png|400px]][[File:FDFD_2D_dielcyl_Im_4.png|400px]]

[[File:FDFD_2D_dielcyl_Re_5.png|400px]][[File:FDFD_2D_dielcyl_Im_5.png|400px]]

[[File:FDFD_2D_benchmark_times_basic.png|800px]]

[[File:FDFD_2D_dx_conv2.png|800px]]

Medusa

2020-08-11T17:36:11Z

BlazS: /* Discussions / Applications */

'''Welcome to the Medusa wiki. To visit the main website, go to [http://e6.ijs.si/medusa/ http://e6.ijs.si/medusa/].'''

In [http://e6.ijs.si/ParallelAndDistributedSystems/ Parallel and Distributed Systems Laboratory] we are working on a C++ library that is first and foremost focused on tools for solving Partial Differential Equations by meshless methods. The basic idea is to create generic codes for tools that are needed for solving not only PDEs but many other problems, e.g. Moving Least Squares approximation, $k$-d tree, domain generation engines, etc.
We call this open source meshless project [http://e6.ijs.si/medusa/ Medusa: Coordinate Free Meshless Method implementation (MM)].

Technical details about code, examples, and can be found on our [http://e6.ijs.si/medusa/docs/ documentation page] and [https://gitlab.com/e62Lab/medusa the code]. [[File:C.png|100px||link=https://gitlab.com/e62Lab/medusa|alt=Alt text|code]] [[File:doxygen.png|100px|link=http://e6.ijs.si/medusa/docs/|alt=Alt text|Documentation page]]

This wiki site is meant for more relaxed discussions about general principles, possible and already implemented applications, preliminary analyses, etc.
Note, that there are many grammatical mistakes, typos, stupid sentences, etc. This wiki is meant for quick information exchange and therefore we do not invest a lot of energy into styling :).

== Documentation ==
* [https://gitlab.com/e62Lab/medusa Code on Gitlab]
* [[How to build | Installation and building]]
* [[Including this library in your project | Including this library in your project]]
* [[Testing | Running tests]]
* [http://e6.ijs.si/medusa/docs/ Technical documentation]
* [[Coding style | Coding style]]
* [[Wiki editing guide | Wiki editing and backup guide]]

== Examples ==
In this section we present exact examples. Each of the below solutions can be found also in in the repository under examples. More explanation about the physical background and solution procedure can be found in following sections.
* [[Philosophy of examples and how to run them]]
* [[Poisson's equation]]
* [[Heat equation]]
* [[Linear elasticity]]
* [[Complex-valued problems]]
* [[Coupled domains]]
* [[Parametric domains | Parametric domains – Curved surface with variable density]]
* [[Customization]]
* [[Ghost nodes]]
* [[Electromagnetic scattering]]
* [[Schrödinger equation]]
* [[Wave equation]]
* [[Cahn-Hilliard equation]]
* [[Meshless Lattice Boltzmann method]]

== Building blocks ==
Medusa is modular coordinate-free parallel implementation of a numerical framework designed, but not limited to, for solving PDEs. In this section we present main modules of the library that can be also used as a standalone tools.
* [[Positioning of computational nodes]]
* [[k-d tree|''k''-d tree]] and other spatial search structures
* Solving [[Solving linear systems | linear systems]], [[Solving overdetermined systems | overdetermined]] and [[Solving underdetermined systems | underdetermined]]
* [[Weighted Least Squares (WLS)]]
* [[Computation of shape functions]]
* [[Meshless Local Strong Form Method (MLSM)]]
* [[Radial basis function-generated finite differences (RBF-FD)]]
* [[Ghost nodes (theory)]]
* [[Integrators for time stepping]]
* [[RBF Interpolation]]

== Discussions / Applications ==
This section is meant for general discussion about the physical background of the examples, the solution procedures, various applications, etc. Note, that code snippets presented in discussion might not reflect the actual state of Medusa.
* Basic PDE solutions
** [[Convection Diffusion equation | Convection Diffusion equation]]
**[[Wave equation application]]
* [[Adaptivity]]
* [[Solid Mechanics]]
** [[Point contact]]
** [[Hertzian contact]]
** [[Cantilever beam]]
** [[Fretting fatigue case]]
* [[Fluid Mechanics]]
** [[Lid driven cavity]]
** [[de Vahl Davis natural convection test]]
** [[Natural convection in 3D irregular domain]]
** [[Natural convection from heated cylinder]]
** [[Natural convection between concentric cylinders]]
* [[Computational electromagnetics]]
** [[Triple dielectric step in 1D]]
** [[Scattering from an infinite cylinder]]
** [[Point source near an anisotropic lens]]
* Other applications
** [[Attenuation due to liquid water content in the atmosphere|Attenuation of a satellite communication]]
** [[Heart rate variability detection]]
** [[Bioheat equation]]

== Performance analyses ==
* [[Execution on Intel® Xeon Phi™ co-processor]]
* [[1D MLSM and FDM comparison]]
* [[:File:tech_report.pdf|Execution overheads due to clumsy types::technical report]] [[File:pdf-file.gif]]
* [[Solving sparse systems]]
* [[Eigen paralelization]]

== Last changes ==
<news unique=1 limit = 5>
*{{{timeanddate}}} :: {{{title}}}

</news>

== Miscellaneous ==
* FAQ - [[Frequently asked questions]].
* [[List of wiki contributors]]
* List of library contributors: [http://e6.ijs.si/medusa/about#about-contributors See the official website]

== References ==
* Slak J., Kosec G. Adaptive radial basis function-generated finite differences method for contact problems. International journal for numerical methods in engineering, ISSN 0029-5981 [http://www-e6.ijs.si/ParallelAndDistributedSystems/pdf/32230439.pdf manuscript]
* Slak J., Kosec G.; Refined meshless local strong form solution of Cauchy-Navier equation on an irregular domain. Engineering analysis with boundary elements. 2018;11 ; [http://comms.ijs.si/~gkosec/data/papers/31107623.pdf manuscript]
* Depolli, M., Kosec, G., Assessment of differential evolution for multi-objective optimization in a natural convection problem solved by a local meshless method. Engineering optimization, 2017, vol. 49, no. 4, pp. 675-692 ;[http://comms.ijs.si/~gkosec/data/papers/29639719.pdf manuscript]
* Kosec G., A local numerical solution of a fluid-flow problem on an irregular domain. Advances in engineering software. 2016;7 ; [29512743] ; [http://comms.ijs.si/~gkosec/data/papers/29512743.pdf manuscript]
* Kosec G., Trobec R., Simulation of semiconductor devices with a local numerical approach. Engineering analysis with boundary elements. 2015;69-75; [27912487] ; [http://comms.ijs.si/~gkosec/data/papers/27912487.pdf manuscript]
* Kosec G., Šarler B., Simulation of macrosegregation with mesosegregates in binary metallic casts by a meshless method. Engineering analysis with boundary elements. 2014;36-44; [http://comms.ijs.si/~gkosec/data/papers/3218939.pdf manuscript]
* Kosec G., Depolli M., Rashkovska A., Trobec R., Super linear speedup in a local parallel meshless solution of thermo-fluid problem. Computers & Structures. 2014;133:30-38; [http://comms.ijs.si/~gkosec/data/papers/27339815.pdf manuscript]
* Kosec G., Zinterhof P., Local strong form meshless method on multiple Graphics Processing Units. Computer modeling in engineering & sciences. 2013;91:377-396; [http://comms.ijs.si/~gkosec/data/papers/26785063.pdf manuscript]
* Kosec G., Šarler B., H-adaptive local radial basis function collocation meshless method. Computers, materials & continua. 2011;26:227-253; [http://comms.ijs.si/~gkosec/data/papers/KosecSarlerBurgers.pdf manuscript]
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321; [http://comms.ijs.si/~gkosec/data/papers/EABE2499.pdf manuscript]
* Kosec G, Zaloznik M, Sarler B, Combeau H. A Meshless Approach Towards Solution of Macrosegregation Phenomena. CMC: Computers, Materials, & Continua. 2011;580:1-27 [http://comms.ijs.si/~gkosec/data/papers/KosecZaloznikSarlerCombeauSegregation.pdf manuscript]
* Kosec G, Sarler B. Solution of thermo-fluid problems by collocation with local pressure correction. International Journal of Numerical Methods for Heat & Fluid Flow. 2008;18:868-82 [http://comms.ijs.si/~gkosec/data/papers/KosecSarlerNS2008.pdf manuscript]
* Trobec R., Kosec G., Parallel Scientific Computing, ISBN: 978-3-319-17072-5 (Print) 978-3-319-17073-2.
* Slak, J., Kosec, G.. Detection of heart rate variability from a wearable differential ECG device., MIPRO 2016, 39th International Convention, 2016, Opatija, Croatia, ISSN 1847-3938, pp 450-455.
* Kolman, M., Kosec, G. Correlation between attenuation of 20 GHz satellite communication link and liquid water content in the atmosphere. MIPRO 2016, 39th International Convention, 2016, Opatija, Croatia, ISSN 1847-3938. pp. 308-313.

==Related pages==
* http://e6.ijs.si/ParallelAndDistributedSystems/#!NumericalMethods
* http://e6.ijs.si/ParallelAndDistributedSystems/#!utils
* http://e6.ijs.si/ParallelAndDistributedSystems/#!NUMA

Computational electromagnetics

2020-08-11T17:34:56Z

BlazS: /* Case studies */

== Case studies ==
The following pages describe the basics of computational electromagnetics, starting with relevant derivations and the basics of classical electromagnetism. The subpages include case studies with analytical solutions if they exist and numerical solutions.
* [[Triple dielectric step in 1D]]
* [[Scattering from an infinite cylinder]]
* [[Point source near an anisotropic lens]]

== Classical electromagnetism ==
=== Maxwell's equations in matter ===
Classical electrodynamics is historically one of the most eminent fields of physics as an extension of classical mechanics, it is very successful in explaining a plethora of phenomena. The dynamics of electric and magnetic fields are described with Maxwell's equations. As we will be studying the interaction of electromagnetic waves with different objects, we need Maxwell's equations in matter
<math>
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\begin{align}
&\nabla \times \b{E}(\b{r}, t) = - \dpar{\b{B}(\b{r}, t)}{t}, \label{eq:TFaraday} \\
&\nabla \times \b{H}(\b{r}, t) = \b{j}(\b{r}, t) + \dpar{\b{D}(\b{r}, t)}{t}, \label{eq:TMaxwell-Ampere} \\
&\nabla \cdot \b{D}(\b{r}, t) = \rho(\b{r}, t), \label{eq:TGaussE} \\
&\nabla \cdot \b{B}(\b{r}, t) = 0. \label{eq:TGaussM}
\end{align}
</math>

The system of equations contains four fields. The electric field strength $\b E$ and density $\b D$ and the magnetic field strength $\b H$ and density $\b B$. The four fields are accompanied by the current density $\b j$ and the charge density $\rho$. For a full description of electromagnetic phenomena, we need to provide another two constitutive equations, that relate the strength and density of the fields
<math>
\newcommand{\eps}{\varepsilon}
\begin{align}
\b{B} &= \mu_0 \mu \b{H}, \label{eq:constM} \\
\b{D} &= \eps_0 \eps \b{E} \label{eq:constE}
\end{align}
</math>
where $\varepsilon_0$ and $\mu_0$ are vacuum permittivity and permeability respectively. The dielectric function $\varepsilon$ and magnetic permeability $\mu$ are in general dependant on both $\b E$ and $\b B$ as well as the frequency $\omega$ and can be second order tensors in anisotropic materials. Equations \eqref{eq:constM} and \eqref{eq:constE} already assume linear material properties, generally the polarisation $\b P$ and magnetisation are defined as power series expansions of the electric field density and magnetic field density, as
<math>
\begin{align}
\label{eq:PMexpans}
\b{P} &= \chi_E \b{D} + \mathcal{O}(D^2), \\
\b{M} &= \chi_M \b{H} + \mathcal{O}(H^2).
\end{align}
</math>
The material linearity assumption holds well for small external fields, meaning small $\b D$ and $\b H$. The treatment of nonlinear terms falls within the field of nonlinear optics and is not relevant for our discussion here.

=== Electromagnetic waves ===
<math>
\def\doubleunderline#1{\underline{\underline{#1}}}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\newcommand{\eps}{\varepsilon}
</math>
When studying light and related phenomena, we are usually interested in the behavior of electromagnetic waves, where the electric and magnetic fields oscillate with constant frequency and propagate in a single direction. The wave formulation helps with understanding optical phenomena such as diffraction and reflection and offers easily graspable quantities such as wave amplitude and wavelength. Wave equations for $\b E$ and $\b H$ can be derived from Maxwell's equations for linear and homogeneous materials. Let's start by treating a more general case, where $\doubleunderline \varepsilon$ is a constant second-order tensor and $\mu$ is a scalar quantity. From Maxwell's equations and constitutive relations we obtain
\begin{equation}
\nabla \times \left[ \doubleunderline{\varepsilon}^{-1} \nabla \times \b{H} \right] = \nabla \times \left[ \doubleunderline{\varepsilon}^{-1} \b j \right] - \mu \mu_0 \varepsilon_0 \ddpar{\b H}{t}.
\end{equation}
and
\begin{equation}
\nabla \times \nabla \times \b{E} = \nabla (\nabla \cdot \b E) - \nabla^2 \b E = -\mu\mu_0\dpar{\b j}{t} - \eps_0 \doubleunderline{\varepsilon} \mu \mu_0 \ddpar{\b E}{t}
\end{equation}
which simplify to wave equations for homogenous (linear) materials and no external currents and charges
\begin{equation}
\nabla^2 \b E = \eps\eps_0\mu\mu_0\ddpar{\b E}{t}, \qquad \text{ or } \qquad \nabla^2 \b H = \eps\eps_0\mu\mu_0\ddpar{\b H}{t}.
\end{equation}
One possible way of solving electromagnetic problems is simply to solve the above wave equations numerically. Another approach is to solve Maxwell's equations in the frequency domain. We start by Fourier expansion of the fields $\b E$, $\b D$, $\b B$ and $\b H$ analogous to
\begin{equation}
\label{eq:Fdecomp}
\b{E}(\b{r}, t) = \int\frac{d\omega}{2\pi} \b{E}(\b{r}, \omega) e^{-i\omega t},
\end{equation}
in the frequency domainall the fields become complex. The ratio between the real and complex component of the fields is the phase difference of the material response to external fields. The Fourier expansion leads to the harmonic form of Maxwell's equations
\begin{equation}
\label{eq:FFaraday}
\nabla \times \b{E}(\b{r}, \omega) = - i\omega \b{B} (\b{r}, \omega),
\end{equation}
\begin{equation}
\label{eq:FMaxwell-Ampere}
\nabla \times \b{H}(\b{r}, \omega) = i\omega \b{D}(\b{r}, \omega).
\end{equation}
The constitutive relations are unchanged. If we assume that there are currents we obtain a wave equation for either $\b D$ or $\b H$,
\begin{equation}
\label{eq:frekaniwave}
\nabla \times \left[ \doubleunderline{\eps}^{-1} \nabla \times \b{H} \right] = \omega^2 \mu_0 \eps_0 \mu \b{H}.
\end{equation}
The above form allows for a calculation of a steady state response to an incident harmonic plane wave, meaning that we can study scattering of any incident wave that can be represented by a plane wave expansion. The time dependent solution is obtained by an inverse Fourier transform. In empty space \eqref{eq:frekaniwave} simplifies to the Helmholtz equation
\begin{equation}
\label{eq:frekemptywave}
\nabla^2 \b{H} = - \omega^2 \mu_0 \eps_0 \b{H} = - k^2 \b{H}.
\end{equation}
Wavenumber $k$ in a vacuum is defined as $k = \omega \sqrt{\mu_0 \epsilon_0} = \frac{2 \pi}{\lambda}$. One of the major advantages of the frequency domain approach over the time domain approach is that the frequency dependence of the dielectric function does not need to be considered, as one simply uses the value of the dielectric function at a certain frequency.

===Boundary conditions for Maxwell's equations===
For modeling purposes the treatment of electromagnetic fields on boundaries between materials with different properties, such as $\varepsilon$ and $\mu$, is crucial. The boundary conditions for Maxwell's equations between two linear and homogenous mediums are as follows
\begin{equation}
\begin{array}{rl}
D_{1 n}-D_{2 n}=\sigma, \qquad & E_{1 t}-E_{2 t}=0 \\
B_{1 n}-B_{2 n}=0, \qquad & H_{1 t}-H_{2 t}=K.\end{array}
\end{equation}
Where we have fields $\b{E}_1, \b{D}_1, \b{B}_1, \b{H}_1$ in medium $1$ and fields $\b{E}_2, \b{D}_2, \b{B}_2, \b{H}_2$ in medium $2$. Subscript $t$ and $n$ denote the tangential and normal component respectively, $\sigma$ is the surface charge density and $\b K$ is the surface current density.

===Frequency dependence of the dielectric function===
As mentioned above the dielectric function $\varepsilon$ is a function of frequency. Some analytical properties of the dielectric function $\varepsilon(\omega)$ are well known, due to their relevance in optics, one example being the Kramers-Kronig relations. The dielectric function is in general, based on a nonlocal response function $\chi$, defined as
\begin{equation}
\varepsilon(\omega) = 1 + \int_{0}^{\infty} \chi(\tau) e^{i \omega \tau} \text{d}\tau.
\end{equation}
For a more detailed analysis of the analytical properties of the function, refer to any standard optics or advanced electromagnetics textbook. When it comes to numerical simulations, models of the dielectric function are of more interest than its analytical properties.
====Dielectric materials====
The simplest model for dielectrics is based on the equation of motion for the bound charges
\begin{equation}
m \frac{\mathrm{d}^{2} \mathbf{r}}{\mathrm{d} t^{2}}=-m \gamma \dot{\mathbf{r}}-m \omega_{0}^{2} \mathbf{r}+e \mathbf{E}(t)
\end{equation}
the terms on the right side describe the dissipation of energy, oscillations and the driving force of the external electric field. A short calculation leads to an expression of the frequency dependence of the dielectric function
\begin{equation}
\varepsilon(\omega)=1+ \frac{\omega_p}{\left(\omega_{0}^{2}-\omega^{2}\right)-i \gamma \omega}.
\end{equation}
The above expression can be simplified for high and low frequencies leading to the so-called Debye and Plasmonic relaxation models. The Lorentz model applies for frequencies where no terms in the denominator can be neglected. Typical Debye and Lorentz relaxations are shown in the plots below.

[[File:Debye.png|400px]][[File:Lorentz.png|400px]]

In general, the dielectric function of a material is modeled as a sum of many Debye and Lorentz relaxations, the parameters of which are determined by microscopic properties of the material itself, an extremely interesting case study is the dielectric function of water.
====Conductors====
When dealing with conductors we use the Drude model of conduction
\begin{equation}
m \frac{d \mathbf{v}(t)}{d t}=-m \gamma \mathbf{v}(t)+e \mathbf{E}(t)
\end{equation}
from which we quickly derive the equation for conductivity
\begin{equation}
\sigma(\omega)=\frac{n e^{2}}{m} \frac{1}{\gamma-i \omega},
\end{equation}
which we can relate to the dielectric function
\begin{equation}
\varepsilon(\omega)=1-\frac{\sigma(\omega)}{i \varepsilon_{0} \omega}.
\end{equation}
==Modeling in the time domain==
The single most frequent type of calculation used in computational electromagnetics is the simulation of electromagnetic fields in the time domain. The calculation propagates $\b E(\b r, t)$ and $\b H(\b r, t)$ in time, usually starting with all fields set to zero and some time dependent current source $\b J(\b r, t)$. The main advantage of time domain approaches is the ability to study transient phenomena while its main drawback is the study resonant problems.
===Finite differences in time domain (FDTD)===
==== Formulation ====
Finite differences in time domain were suggested in 1966 by Yee. The method uses finite differences to discretize Maxwell's equations and propagates the fields in time by alternating between updating the electric field strength and magnetic field strength. One of its defining attributes is the use of a staggered grid to represent $\b H$ and $\b E$ nodes. The staggered grid simplifies the calculation of curl expressions and helps with numerical stability and treating the boundaries of materials.

[[File:FDTD2.png|400px|center]]

This type of solution procedure provides an alternative to solving the wave equation for either $\b H$ or $\b E$.

==== Modeling sources ====
Since we are starting with all fields set to zero, we need to somehow introduce power into the model. This is usually done with current sources. The current signal used depends on the type of calculation. If we are interested in a steady-state solution we would use a sinusoidal source, that is slowly turned on and simulated the device until a steady state is reached. On the other hand, if we are interested in calculating the response of a device to many frequencies, we would use a Gaussian pulse and simulate the device until all the energy has left the system.

Taflove in his book [1] differentiates between five types of sources. The easiest to implement are so-called ''hard sources'', where one simply sets the value of a field in a node to the source value e.g.
\begin{equation}
\b E_z(i, t) = \sin 2\pi ft.
\end{equation}
The problem is that the node now acts as a perfect dielectric and incoming waves reflect from the point and do not pass through as they should. On the other hand one can use a ''soft source'' where energy is inserted into the domain via electric and magnetic (non-physical) current
\begin{equation}
\b E_z(i, t) = E_z(i, t) + J(t) \qquad \text{in} \qquad \b H_x(i, t) = H_x(i, t) + M(t).
\end{equation}
Waves no longer reflect from the source, however, the source still radiates in all directions, which is sometimes not desired. This issue is solved with TFSF (total field, scattered field) sources, which are based on the linearity of Maxwell's equations and allow to construct sources that inject the energy only in the desired direction. This is accomplished by decomposing the total fields into the scattered and incident
\begin{equation}
\b{E} = \b{E}^i + \b{E}^s \text{, } \qquad \b{H} = \b{H}^i + \b{H}^s.
\end{equation}
and then dividing the computational domain $\Omega$ into two subdomains. In one of the subdomains, we calculate only using the scattered part of the fields and in the other, we calculate with the total field. With careful treatment of the boundary between the subdomains one can construct an arbitrary incident wave. Alongside the already mentioned source types, ''waveguide'' sources and ''total scattered field'' sources. For an in-depth discussion refer to Taflove [1].

==== Transmittance and Reflectivity ====
Transmittance and reflectivity are dependent on the frequency of the incoming wave. One of the main strengths of time domain calculations is that one can calculate both coefficients for a wide frequency range with a single simulation. This is achieved by using a source that contains a wide range of frequencies i.e. a short Gaussian pulse.

Let us illustrate how to calculate $R$ and $T$ in one dimension. The calculation is similar in higher dimensions. We make two assumptions, firstly the source is directed directly into the device. This means that we will have only the reflected field on one side of the device and only the transmitted field on the other side of the device. Secondly, we assume that nothing reflects from the edges of the domain. We record the signals on both edges of the domain and label them $E_r(t)$ for reflected and $E_t(t)$ for transmitted. It is important to simulate for a long enough time that there is no energy left in the device. Otherwise, the calculation of the coefficients has not converged yet. The Fourier transforms of the signals at the edges, along with the transform of the source, define $R$ and $T$
\begin{equation}
R(\omega) = \left| \frac{E_o(\omega)}{S(\omega)} \right| \qquad \text{and} \qquad T(\omega)= \left| \frac{E_p(\omega)}{S(\omega)} \right|.
\end{equation}
Since we only know the field values at discrete time intervals the transformation is discrete
\begin{equation}
F(\omega) \cong \Delta t \sum_{m=1}^{M}\left(e^{-j \omega \Delta t}\right)^{m} \cdot f(m \Delta t),
\end{equation}
in practice, so that it is not needed to store the signals on both edges which may take up considerable computer memory, the sum can be calculated after each step in the simulation to save space.

<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}
</math>

==Modeling in the frequency domain==
In the frequency domain, two types of calculations are most common. Within the field of photonic crystals, bandstructures of periodic materials are sought by solving a generalized eigenvalue problem $Ax = \omega^2Bx$. In engineering applications, a response of a structure to a periodical source $\b J(x) e^{-i\omega t}$ is more often calculated. We will focus on the second type of calculation.

==== Formulation ====
We need to solve the equation that can be derived by combining Maxwell's equations and constitutive relations
\begin{equation}
\left[(\nabla \times \nabla \times)-\frac{\omega^{2}}{c^{2}} \varepsilon(\mathbf{r})\right] \mathbf{E}(\mathbf{r})=i \omega \mu_{0} \mathbf{J}(\mathbf{r}).
\end{equation}
The electric field strength $\b E(\b r)$ is complex, the time dependent solution is obtained via $\b E(\b r,t) = \b E(\b r) e^{i\omega t}$. The magnetic field can be obtained by taking a curl of the solution above
\begin{equation}
\mathbf{H}=-\frac{i}{\omega \mu_{0}} \nabla \times \mathbf{E}.
\end{equation}

==== Absorbing boundary conditions ====
Absorbing boundary conditions in the frequency domain can be realized in a number of ways. One of the most straight forward approaches is to impose the Sommerfeld radiation condition or a numerical approximation of it. However, it is more common to surround the computational domain with a thin absorbing layer, which does not reflect outgoing waves. The perfectly matched layer (PML) is one of such absorbing boundary conditions. There are many variants of PML's, but it is clear that the use of one worsens the conditioning of the problem matrix [2]. As it turns out UPML (Uniaxial PML) is numerically less stable than SC-PML (stretched coordinate PML), but efficient preconditioners exist for cases when SC-PML is hard to implement.

===== Uniaxial Perfectly Matched Layer =====
Inside the computational domain we have an amplitude equation for $\b E$,
\begin{equation}
\nabla \times \doubleunderline{\mu}^{-1} \nabla \times \mathbf{E}-\omega^{2} \varepsilon \mathbf{E}=-\mathrm{i} \omega \mathbf{J}.
\end{equation}
We surround the computational domain $\Omega$ with a thin domain $\Omega'$, where we have an uniaxial anisotropic material, which damps the electromagnetic field
\begin{equation}
\nabla \times \doubleunderline{\mu_{s}}^{-1} \nabla \times \mathbf{E}-\omega^{2} \doubleunderline{\varepsilon_{s}} \mathbf{E}=-\mathrm{i} \omega \mathbf{J}
\end{equation}
dielectric and permeability tensors are diagonal,
\begin{equation}
\doubleunderline{\varepsilon_{s}}=\varepsilon\left[\begin{array}{ccc}\frac{s_{y} s_{z}}{s_{x}} & 0 & 0 \\ 0 & \frac{s_{s} s_{x}}{s_{y}} & 0 \\ 0 & 0 & \frac{s_{x} s_{y}}{s_{z}}\end{array}\right], \quad \doubleunderline{\mu_{s}}=\mu\left[\begin{array}{ccc}\frac{s_{y} s_{z}}{s_{x}} & 0 & 0 \\ 0 & \frac{s_{z} s_{x}}{s_{y}} & 0 \\ 0 & 0 & \frac{s_{x} s_{y}}{s_{z}}\end{array}\right].
\end{equation}
To avoid reflections of waves from the boundary between layers, we slowly turn on the loss parameter $s_w$
\begin{equation}
s_{w}(l)=1-\mathrm{i} \frac{\sigma_{w}(l)}{\omega \varepsilon_{0}}
\end{equation}
where $l$ is the distance from the boundary of the PML.

===== Streched Coordinate Perfectly Matched Layer =====
In the stretched coordinate PML the coordinates are stretched and the PML no longer resembles any physical material. The equation inside the PML remains unchanged,
\begin{equation}
\label{eq:SC-PML}
\nabla_{s} \times \mu^{-1} \nabla_{s} \times \mathbf{E}-\omega^{2} \varepsilon \mathbf{E}=-\mathrm{i} \omega \mathbf{J}.
\end{equation}
However, the gradient is redefined as
\begin{equation}
\nabla_{s}=\hat{e}_x \frac{1}{s_{x}} \frac{\partial}{\partial x}+\hat{e}_y \frac{1}{s_{y}} \frac{\partial}{\partial y}+\hat{e}_z \frac{1}{s_{z}} \frac{\partial}{\partial z}.
\end{equation}
Just like with UPML, the damping is slowly turned on. SC-PML is considered more stable than UMPL, but it is usually trickier to implement. For a magnetically homogenous and isotropic medium the SC-PML equation simplifies to
\begin{equation}
-\nabla^2_s\b E + \nabla_s\left( \nabla_{s} \cdot \b E \right) -\omega^2 \eps \mu \b E = -i\omega\b J,
\end{equation}
in matrix form
\begin{multline}
-
\begin{bmatrix}
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_x}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_x}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_x}{z}\right) \\
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_y}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_y}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_y}{z}\right) \\
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_z}{x}\right) +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_z}{y}\right) +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_z}{z}\right)
\end{bmatrix}
+
\begin{bmatrix}
\frac{1}{s_x} \frac{\partial}{\partial x} \left(\frac{1}{s_x}\dpar{E_x}{x}\right) +
\frac{1}{s_y s_x}\frac{\partial E_y}{\partial x \partial y} +
\frac{1}{s_z s_x}\frac{\partial E_z}{\partial z \partial z} \\
\frac{1}{s_x s_y}\frac{\partial E_x}{\partial y \partial x} +
\frac{1}{s_y} \frac{\partial}{\partial y} \left(\frac{1}{s_y}\dpar{E_y}{y}\right) +
\frac{1}{s_z s_y}\frac{\partial E_z}{\partial y \partial z} \\
\frac{1}{s_x s_z}\frac{\partial E_x}{\partial z \partial x} +
\frac{1}{s_y s_z}\frac{\partial E_y}{\partial z \partial y} +
\frac{1}{s_z} \frac{\partial}{\partial z} \left(\frac{1}{s_z}\dpar{E_z}{z}\right)
\end{bmatrix} -
\\ -\omega^2 \eps \mu \b E
\begin{bmatrix}
E_x \\
E_y \\
E_z
\end{bmatrix}
= -i\omega
\begin{bmatrix}
J_x \\
J_y \\
J_z
\end{bmatrix}
\end{multline}
The expression further simplifies in cases where some of $s_w$'s are constant (cartesian coordinates).
===== Benchmarking =====
Parameters of the PML are usually determined empirically, separately for each simulation, however, there are some rules of thumb. PMLs by Berenger's definition are exact for a continuous equation. When the equation is discretized reflections appear at the boundary of the PML. A lot of care needs to be taken so that the reflections are minimized. This is usually achieved by giving $\sigma_{w}(l)$ a polynomial relationship with the distance $l$ from the PML boundary
\begin{equation}
\sigma_{w}(l,) = \sigma_{max}\left(\frac{l}{d}\right)^m
\end{equation}

where $\sigma_{max}$ is the maximum loss parameter, $d$ is the thickness of the layer and $m$ is the power of the polynomial, and usually lies somewhere between $2$ and $4$. For a normally incident wave on the PML, we can choose the Reflectance of the PML and accordingly tune parameters $d$ and $m$ via
\begin{equation}
\sigma_{w, \max }=-\frac{(m+1) \ln R}{2 \eta_{0} d}
\end{equation}
where $\eta_0$ is the Impedance of free space.

==References==
# Taflove, A. and Hagness, S.C., 2005. Computational electrodynamics: the finite-difference time-domain method. Artech house.
# Shin, W. and Fan, S., 2012. Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell’s equations solvers. Journal of Computational Physics, 231(8), pp.3406-3431.
# Podgornik, Rudolf, and Andrej Vilfan. Elektromagnetno polje. DMFA-založništvo, 2012.
# John, S. G. J., et al. "Photonic crystals: molding the flow of light." Princeton University of press. 2008.

Scattering from an infinite cylinder

2020-08-11T17:34:37Z

BlazS: /* Results */

Scattering from an infinite cylinder

2020-08-11T17:33:24Z

BlazS: /* Case */

Scattering from an infinite cylinder

2020-08-11T17:32:48Z

BlazS: /* Case */

Scattering from an infinite cylinder

2020-08-11T17:31:46Z

BlazS: /* Case */

Back to [[Computational electromagnetics]].

==Case==
As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. Zgled je pomemben iz dveh razlogov, prvi\v c, ogledali si bomo kako vpadno polje zapi\v semo v rob med materialom in prostim prostorom in lahko kot take vire analiti\v cno predstavimo\footnote{Tu to ne pride posebno prav, ker imamo opravka z ravnim valom, v primerih ko imamo to\v ckaste izvire pa je zelo prakti\v cno.}. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
\begin{align}
\nabla^2 v +\eps_r \frac{\omega^2}{c_0^2}\thinspace v = 0 \qquad &\textnormal{in} \quad \Omega \label{eq:inner} \\
\nabla^2 u^s + \frac{\omega^2}{c_0^2}u^s = 0 \qquad & \textnormal{outside} \quad \Omega
\label{eq:outer}
\end{align}
The field $v$ is the total field inside, $u^s$ is only the scattered part of $E_z$ Outside the cylinder. The boundary conditions are
\begin{equation}
v - u^s =u^i \qquad \textnormal{on} \quad \partial \Omega \label{eq:BC1}, \qquad \text{in} \qquad
\dpar{v}{n} - \frac{1}{\eps_r}\dpar{u^s}{n} = \dpar{u^i}{n} \qquad \textnormal{on} \quad \partial \Omega.
\end{equation}
The incident field $u^i$ has an analytical form in our case a plane wave $e^{ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.

==Analytical solution==
==Results==
[[File:FDFD_2D_dielcyl_Im_1.png|400px]]
[[File:FDFD_2D_dielcyl_Im_2.png|400px]]
[[File:FDFD_2D_dielcyl_Im_3.png|400px]]
[[File:FDFD_2D_dielcyl_Im_4.png|400px]]
[[File:FDFD_2D_dielcyl_Im_5.png|400px]]
[[File:FDFD_2D_dielcyl_Re_1.png|400px]]
[[File:FDFD_2D_dielcyl_Re_2.png|400px]]
[[File:FDFD_2D_dielcyl_Re_3.png|400px]]
[[File:FDFD_2D_dielcyl_Re_4.png|400px]]
[[File:FDFD_2D_dielcyl_Re_5.png|400px]]
[[File:FDFD_2D_benchmark_times_basic.png|400px]]
[[File:FDFD_2D_dx_conv2.png|400px]]

Scattering from an infinite cylinder

2020-08-11T17:20:10Z

BlazS:

Back to [[Computational electromagnetics]].

==Case==
==Analytical solution==
==Results==
[[File:FDFD_2D_dielcyl_Im_1.png|400px]]
[[File:FDFD_2D_dielcyl_Im_2.png|400px]]
[[File:FDFD_2D_dielcyl_Im_3.png|400px]]
[[File:FDFD_2D_dielcyl_Im_4.png|400px]]
[[File:FDFD_2D_dielcyl_Im_5.png|400px]]
[[File:FDFD_2D_dielcyl_Re_1.png|400px]]
[[File:FDFD_2D_dielcyl_Re_2.png|400px]]
[[File:FDFD_2D_dielcyl_Re_3.png|400px]]
[[File:FDFD_2D_dielcyl_Re_4.png|400px]]
[[File:FDFD_2D_dielcyl_Re_5.png|400px]]
[[File:FDFD_2D_benchmark_times_basic.png|400px]]
[[File:FDFD_2D_dx_conv2.png|400px]]

File:FDFD 2D dx conv2.png

2020-08-11T17:19:03Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

File:FDFD 2D benchmark times basic.png

2020-08-11T17:19:02Z

BlazS: File uploaded with MsUpload

File uploaded with MsUpload

Point source near an anisotropic lens

2020-08-11T17:02:31Z

BlazS: /* Results */

File:FDFD 2D anizolens Im.png

2020-08-11T17:02:03Z

BlazS: BlazS uploaded a new version of File:FDFD 2D anizolens Im.png

File uploaded with MsUpload

File:FDFD 2D anizolens Re.png

2020-08-11T17:02:03Z

BlazS: BlazS uploaded a new version of File:FDFD 2D anizolens Re.png

File uploaded with MsUpload

Point source near an anisotropic lens

2020-08-11T17:00:48Z

BlazS: /* Case */

Point source near an anisotropic lens

2020-08-11T16:11:51Z

BlazS:

Back to [[Computational electromagnetics]].

==Case==
==Results==
<gallery>
File:FDFD_2D_anizolens_Im.png
File:FDFD_2D_anizolens_Re.png
</gallery>