Electromagnetics Mini Applications

$\newcommand{\A}{\vec{A}}\newcommand{\B}{\vec{B}} \newcommand{\D}{\vec{D}}\newcommand{\E}{\vec{E}} \newcommand{\H}{\vec{H}}\newcommand{\J}{\vec{J}} \newcommand{\M}{\vec{M}}\newcommand{\P}{\vec{P}} \newcommand{\F}{\vec{F}} \newcommand{\dd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\cross}{\times}\newcommand{\inner}{\cdot} \newcommand{\div}{\nabla\cdot}\newcommand{\curl}{\nabla\times} \newcommand{\grad}{\nabla}$

The miniapps/electromagnetics directory contains a collection of electromagnetic miniapps based on MFEM.

Compared to the example codes, the miniapps are more complex, demonstrating more advanced usage of the library. They are intended to be more representative of MFEM-based application codes. We recommend that new users start with the example codes before moving to the miniapps.

The current electromagnetic miniapps are described below.

Electromagnetics

The equations describing electromagnetic phenomena are known collectively as the Maxwell Equations. They are usually given as:

$$\begin{align} \curl\H - \dd{\D}{t} & = \J \label{ampere} \\ \curl\E + \dd{\B}{t} & = 0 \label{faraday} \\ \div\D & = \rho \label{gauss} \\ \div\B & = 0 \label{divb} \end{align}$$

Where equation \eqref{ampere} can be referred to as Ampère's Law, equation \eqref{faraday} is called Faraday's Law, equation \eqref{gauss} is Gauss's Law, and equation \eqref{divb} doesn't generally have a name but is related to the nonexistence of magnetic monopoles. The various fields in these equations are:

Symbol	Name	SI Units
$\H$	magnetic field	Ampere/meter
$\B$	magnetic flux density	Tesla
$\E$	electric field	Volt/meter
$\D$	electric displacement	Coulomb/meter$^2$
$\J$	current density	Ampere/meter$^2$
$\rho$	charge density	Coulomb/meter$^3$

In the literature these names do vary, particularly those for $\H$ and $\B$, but in this document we will try to adhere to the convention laid out above.

Generally we also need constitutive relations between $\E$ and $\D$ and/or between $\H$ and $\B$. These relations start with the definitions:

$$\begin{align} \D & = \epsilon_0\E + \P \label{const_d} \\ \B & = \mu_0(\H + \M) \label{const_b} \end{align}$$

Where $\P$ is the polarization density, and $\M$ is the magnetization. Also, $\epsilon_0$ is the permittivity of free space and $\mu_0$ is the permeability of free space which are both constants of nature. In many common materials the polarization density can be approximated as a scalar multiple of the electric field, i.e., $\P = \epsilon_0\chi\E$, where $\chi$ is called the electric susceptibility. In such cases we usually use the relation $\D = \epsilon\E$ with $\epsilon = \epsilon_0(1 + \chi)$ and call $\epsilon$ the permittivity of the material.

The nature of magnetization is more complicated but we will take a very simplified view which is valid in many situations. Specifically, we will assume that either $\M$ is proportional to $\H$ yielding the relation $\B = \mu\H$ where $\mu = \mu_0(1 + \chi_M)$ and $\chi_M$ is the magnetic susceptibility or that $\M$ is independent of the applied field. The former case pertains to both diamagnetic and paramagnetic materials and the latter to ferromagnetic materials.

Finally we should note that equations \eqref{ampere} and \eqref{gauss} can be combined to yield the equation of charge continuity $\dd{\rho}{t} + \div\J = 0$ which can be important in plasma physics and magnetohydrodynamics (MHD).

Electrostatics

Electrostatic problems come in a variety of subtypes but they all derive from Gauss's Law and Faraday's Law (equations \eqref{gauss} and \eqref{faraday}). When we assume no time variation, Faraday's Law becomes simply $\curl\E = 0$. This suggests that the electric field can be expressed as the gradient of a scalar field which is traditionally taken to be $-\varphi$, i.e.

$$\E = -\grad\varphi \label{gradphi}$$

where $\varphi$ is called the electric potential and has units of Volts in the SI system. Inserting this definition into equation \eqref{gauss} gives:

$$-\div\epsilon\grad\varphi = \rho - \div\P \label{poisson}$$

which is Poisson's equation for the electric potential, where we have assumed a linear constitutive relation between $\D$ and $\E$ of the form $\D = \epsilon\E + \P$. This allows a polarization which is proportional to $\E$ as well as a polarization independent of $\E$. If this relation happens to be nonlinear then Poisson's equation would need to be replaced with a more complicated nonlinear expression.

The solutions to equation \eqref{poisson} are non unique because they can be shifted by any additive constant. This means that we must apply a Dirichlet boundary condition at least at one point in the problem domain in order to obtain a solution. Typically this point will be on the boundary but it need not be so. Such a Dirichlet value is equivalent to fixing the voltage (a.k.a. potential) at one or more locations. Additionally, this equation admits a normal derivative boundary condition. This corresponds to setting $\hat{n}\cdot\D$ to a prescribed value on some portion of the boundary. This is equivalent to defining a surface charge density on that portion of the boundary.

Volta Mini Application

The electrostatics mini application, named volta after the inventor of the voltaic pile, is intended to demonstrate how to solve standard electrostatics problems in MFEM. Its source terms and boundary conditions are simple but they should indicate how more specialized sources or boundary conditions could be implemented.

Note that this application assumes the mesh coordinates are given in meters.

Mini Application Features

Permittivity: The permittivity, $\epsilon$, is assumed to be that of free space except for an optional sphere of dielectric material which can be defined by the user. The command line option -ds can be used to set the parameters for this dielectric sphere. For example, to produce a sphere at the origin with a radius of 0.5 and a relative permittivity of 3 the user would specify: -ds '0 0 0 0.5 3'.

Charge Density: The charge density, $\rho$, is assumed to be zero except for an optional sphere of uniform charge density which can be defined by the user. The command line option for this is -cs which follows the same pattern as the dielectric sphere. Note that the last entry is the total charge of the sphere and not its charge density.

Polarization: A polarization vector function, $\P$, can be imposed as a source of the electric field. The command line option -vp creates a polarization due to a simple voltaic pile, i.e., a cylinder which is electrically polarized along its axis. The user should specify the two end points of the cylinder axis, its radius and the magnitude of the polarization vector.

Dirichlet BC: Dirichlet boundary conditions can either specify piecewise constant voltages on a collection of surfaces or they can specify a gradient field which approximates a uniform applied electric field. In either case the user specifies the surfaces where the Dirichlet boundary condition should be applied using the -dbcs option followed by a list of boundary attributes. For example to select surfaces 2, 3, and 4 the user would use the following: -dbcs '2 3 4'.

To apply a gradient field on these surfaces the user would also use the -dbcg option. This defaults to the uniform field $\E = (0,0,1)$ in 3D or $\E = (0,1)$ in 2D. An arbitrary vector can be specified with -uebc followed by the desired vector, e.g., to apply $\E = (1,2,3)$ the user would supply: -uebc '1 2 3'.

To specify piecewise constant potential values the user would list the desired values after -dbcv as follows: -dbcv '0.0 1.0 -1.0'.

Neumann BC: Neumann boundary conditions set the normal component of the electric displacement on portions of the boundary. This normal component is equivalent to the surface charge density on the surface. This is rarely used because surface charge densities are rarely known unless they are known to be zero. However, if the surface charge density is zero then the Neumann BCs are not needed because this is the natural boundary condition. Only piecewise constant Neumann boundary conditions are supported. They can be set analogously to piecewise Dirichlet boundary conditions but using options -nbcs and -nbcv.

Magnetostatics

Magnetostatic problems arise when we assume no time variation in Ampère's Law \eqref{ampere} which leads to:

$$\curl\H = \J \nonumber$$

We will again assume a somewhat more general constitutive relation between $\H$ and $\vec{B}$ than is normally seen:

$$\B = \mu\H + \mu_0\M = \mu_0(1 + \chi_M)\H + \mu_0\M \nonumber$$

Where the magnetization is split into two portions; one which is proportional to $\H$ and given by $\chi_M\H$, and another which is independent of $\H$ and is given by $\M$. This allows for paramagnetic and/or diamagnetic materials defined through $\mu$ as well as ferromagnetic materials represented by $\M$. This choice yields:

$$\curl\mu^{-1}\B = \J + \curl\mu^{-1}\mu_0\M \nonumber$$

Which, when combined with equation \eqref{divb}, becomes:

$$\curl\mu^{-1}\curl\A = \J + \curl\mu^{-1}\mu_0\M $$

If $\J$ happens to be zero we have another option because we can assume that $\H = -\grad\varphi_M$ for some scalar potential $\varphi_M$. When combined with equation \eqref{divb} this leads to:

$$\div\mu\grad\varphi_M = \div\mu_0\M $$

Currently only the vector potential equation is used so we will focus on that for the remainder of this document.

The vector potential is again non unique so we must apply additional constraints in order to arrive at a solution for $\A$. When working analytically it is common to constrain the solution by restricting the divergence of $\A$ but numerically this leads to other complications. For our problems of interest it will be necessary to require Dirichlet boundary conditions on the entire outer surface in order to sufficiently constrain the solution.

Dirichlet boundary conditions for the vector potential on a surface provide a means to specify the component of $\B$ normal to that surface. For example, setting the tangential components of $\A$ to be zero on a particular surface results in a magnetic flux density which must be tangent to that surface.

Tesla Mini Application

The magnetostatics mini application, named tesla after the unit of magnetic field strength (and of course the man Nikola Tesla), is intended to demonstrate how to solve standard magnetostatics problems in MFEM. Its source terms and boundary conditions are simple but they should indicate how more specialized sources of boundary conditions could be implemented. A detailed overview of the equations being solved and their discretization can be found here: Tesla Theory Notes.

Note that this application assumes the mesh coordinates are given in meters.

Mini Application Features

Permeability: The permeability, $\mu$, is assumed to be that of free space except for an optional spherical shell of diamagnetic or paramagnetic material which can be defined by the user. The command line option -ms can be used to set the parameters for this shell.

For example, to produce a shell at the origin with inner and outer radii of 0.4 and 0.5 respectively and a relative permeability of 3 the user would specify: -ms '0 0 0 0.4 0.5 3'.

Current Density: The current density, $\J$, is assumed to be zero except for an optional ring of constant current which can be defined by the user. The command line option for this is -cr which requires two points giving the end points of the ring's axis, inner and outer radii, and a constant total current.

For example, to specify a ring centered at the origin and laying in the XY plane with a thickness of 0.2 and radii 0.4 and 0.5, and a current of 2 amps the user would give: -cr 0 0 -0.1 0 0 0.1 0.4 0.5 2.

Magnetization: A permanent magnetization, $\M$, can be applied in the form of a cylindrical magnet with poles at its circular ends. The command line option is -bm which indicates a 'bar magnet'. The option requires the two end points of the cylinder's axis, its radius, and the magnitude of the magnetization.

Surface Current Density: A surface current can be imposed indirectly by specifying separate surface patches with different voltages as well as a collection of surface patches connecting the voltages through which the current will flow. The voltage surfaces and their voltages can be specified using -vbcs followed by the indices of the surfaces and -vbcv followed by their voltages. The path for the surface current ($\vec{K}$) is specified by using -kbcs followed by a set of surface indices.

For example, applying voltages 1 and -1 to surfaces 2 and 3 with a current path along surfaces 4 and 6 would be specified as: -vbcs '2 3' -vbcv '1 -1' -kbcs '4 6'.

Any surfaces not listed as voltage or current surfaces will be assigned as homogeneous Dirichlet boundaries. Note that when this option is selected an auxiliary electrostatic problem will be solved on the surface of the geometry to compute the surface current.

Dirichlet BC: Dirichlet boundary conditions are required if a surface current density is not defined. For this reason the user need not specify boundary surfaces by number since the boundary condition must be applied on all of them. The default boundary condition is a homogeneous Dirichlet boundary condition on all outer surfaces. This means that the normal component of $\B$ will be zero at the outer boundary. An alternative is to specify a desired uniform magnetic flux density on the entire outer surface. This is accomplished with the -ubbc command line option followed by the desired $\B$ vector.

Transient Full-Wave Electromagnetics

Transient electromagnetics problems are governed by the time-dependent Maxwell equations \eqref{ampere} and \eqref{faraday} when combined using the constitutive relations \eqref{const_d} and \eqref{const_b}. When combined these equations can describe the evolution and propagation of electromagnetic waves.

$$\begin{align} \dd{(\epsilon\E)}{t} & = \curl(\mu^{-1}\B) - \sigma \E - \J \\ \dd{\B}{t} & = - \curl\E \end{align}$$

The term $\sigma\E$ arises in the presence of electrically conductive materials where the electric field induces a current which can be separated from $\J$. In such cases the total current appearing in Ampère's Law \eqref{ampere} can be expressed as the sum of an applied current (also labeled as $\J$) and an induced current $\sigma\E$.

Solving these equations requires initial conditions for both the electric and magnetic fields $\E$ and $\B$ as well as boundary conditions related to the tangential components of $\E$ or $\H$. Other formulations are possible such as evolving $\H$ and $\D$ or the potentials $\varphi$ and $\A$. This system of equations can also be written as a single second order equation involving only $\E$, $\H$, $\varphi$, or $\A$. Each of these formulations has a different set of sources, initial and boundary conditions for which it is well-suited. The choice we make here is perhaps the most common but it may not be the most convenient choice for a given application.

These equations can be used to evolve their initial conditions or they can be driven by either a current source or through time-varying boundary conditions. It is also possible to combine all three of these sources in a single simulation.

Maxwell Mini Application

The electrodynamics mini application, named maxwell after James Clerk Maxwell who first formulated the classical theory of electromagnetic radiation, is intended to demonstrate how to solve transient wave problems in MFEM. Its source terms and boundary conditions are simple but they should indicate how more specialized sources or boundary conditions could be implemented. A detailed overview of the equations being solved and their discretization can be found here: Maxwell Theory Notes.

An example simulation is depicted below (click to animate the wave propagation).

Time integration is handled by a variable order symplectic time integration algorithm. This algorithm is designed for systems of equations which are derived from a Hamiltonian and it helps to ensure energy conservation within some tolerance. The time step used during integration is automatically chosen based on the largest stable time step as computed from the largest eigenvalue of the update equations. This determination involves a user-adjustable factor which creates a safety margin. By default the actual time step is less than 95% of the estimate for the largest stable time step.

Note that this application assumes the mesh coordinates are given in meters. Internally the code assumes time is in seconds but the command line options use nanoseconds for convenience.

Mini Application Features

Time Evolution: The initial and final times for the simulation can be specified, in nanoseconds, with the -ti and -tf options. Visualization snapshots of data will be written out after time intervals specified by -ts which again given in nanoseconds. The order of the time integration can be specified, from 1 to 4, using the -to option.

For example, to produce a shell at the origin with inner and outer radii of 0.4 and 0.5 respectively and a relative permeability of 3 the user would specify: -ms '0 0 0 0.4 0.5 3'.

Conductivity: The conductivity, $\sigma$, is assumed to be zero except for an optional sphere of conductive material which can be defined by the user. The command line option -cs can be used to set the parameters for this conductive sphere. For example, to produce a sphere at the origin with a radius of 0.5 and a conductivity of 3,000,000 S/m the user would specify: -cs '0 0 0 0.5 3e6'.

Current Density: The current density, $\J$, is assumed to be zero except for an optional cylinder of pulsed current which can be defined by the user. The command line option for this is -dp, short for 'dipole pulse', which requires two points giving the end points of the cylinder's axis, radius, amplitude ($\alpha$), pulse center ($\beta$), and a pulse width ($\gamma$). The time dependence of this pulse is given by: $$\J(t) = \hat{a} \alpha e^{-(t-\beta)^2/(2\gamma^2)}$$ Where $\hat{a}$ is the unit vector along the cylinder's axis and both $\beta$ and $\gamma$ are specified in nanoseconds.

Dirichlet BC: Homogeneous Dirichlet boundary conditions, which constrain the tangential components of $\frac{\partial\E}{\partial t}$ to be zero, can be activated on a portion of the boundary by specifying a list of boundary attributes such as -dbcs '4 8'. For convenience a boundary attribute of '-1' can be used to specify all boundary surfaces. Non-Homogeneous, time-dependent Dirichlet boundary conditions are supported by the Maxwell solver so a user can edit maxwell.cpp and supply their own function if desired.

Absorbing BC: A first order Sommerfeld absorbing boundary condition can be applied to a portion of the boundary using the -abcs option along with a list of boundary attributes such as -abcs '4 18'. Again, the special purpose boundary attribute '-1' can be used to specify all boundary surfaces. This boundary condition depends on a coefficient, $\eta^{-1}=\sqrt{\epsilon/\mu}$, which must be matched to the materials just inside the boundary. The code assumes that the permittivity and permeability are those of the vacuum near the surface but, if this is not the case, an ambitious user can replace etaInvCoef_ with a more appropriate function.

Transient Magnetics and Joule Heating

Joule Mini Application

The transient magnetics mini application, named joule after the SI unit of energy (and the scientist James Prescott Joule, who was also a brewer), is intended to demonstrate how to solve transient implicit diffusion problems. The equations of low-frequency electromagnetics are coupled with the equations of heat transfer. The coupling is one way, electromagnetics generates Joule heating, but the heating does not affect the electromagnetics. The thermal problem is solved using an $H(\mathrm{div})$ method, i.e. temperature is discontinuous and the thermal flux $\F$ is in $H(\mathrm{div})$. There are three linear solves per time step:

Poisson's equation for the scalar electric potential is solved using the AMG preconditioner,
the electric diffusion equation is solved using the AMS preconditioner, and
the thermal diffusion equation is solved using the ADS preconditioner.

Two example meshes are provided, one is a straight circular metal rod in vacuum, the other is a helical coil in vacuum (the latter is 21MB and can be downloaded from here). The idea is that a voltage is applied to the ends of the rod/coil, the electric field diffuses into the metal, the metal is heated by Joule heating, the heat diffuses out.

The equations are:

$$\begin{align} \div\sigma\grad\Phi &= 0 \\ \sigma \E &= \curl\mu^{-1} \B - \sigma \grad \Phi \\ \frac{d \B}{d t} &= - \curl \E \\ \F &= -k \grad T \\ c \frac{d T}{d t} &= - \div \F + \sigma \E \cdot \E \end{align}$$

The equations are integrated in time using implicit time integration, either midpoint or higher order SDIRK.

Since there are three solves, three sets of boundary conditions must be specified. The essential BC's are the scalar potential, the electric field, and the thermal flux. These are not set via command line arguments, you have to edit the code to change these. To change these, search the code for ess_bdr.

There are conducting and non-conducting material regions, and the mesh must have integer attributes to specify these regions. To change these, search the code for std::map<int, double> this maps the integer attribute to the floating-point material value.

Note that this application assumes the mesh coordinates are given in meters.

The above picture shows Joule heating of a cylinder using the mesh cylinder-hex.mesh. The cylinder is surrounded by vacuum. The black arrows show the magnetic field $\B$, the magenta arrows show the heat flux $\F$, and the pseudocolor in the center of the cylinder shows the temperature.

Mini Application Features

Boundary Conditions: Since there are three solves, three sets of boundary conditions must be specified. The essential BC's are the voltage for the scalar potential, the tangential electric field, and the normal thermal flux. These are not set via command line arguments, you have to edit the code to change these. To change these, search the code for ess_bdr. Note that the essential BC's can be time varying.

Material Properties: There are conducting and non-conducting material regions, and the mesh must have integer attributes to specify these regions. To change these, search the code for std::map<int, double> this maps the integer attribute to the floating-point material value.