 |
deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00:01+00:00
|
Loading...
Searching...
No Matches
|
deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00:01+00:00
|
This tutorial depends on step-97.
| Table of contents | |
|---|---|
| |
This program was contributed by Siarhei Uzunbajakau, CEM Books, 2026.
This tutorial illustrates the application of the FE_Q, FE_Nedelec, FE_RaviartThomas, and FE_DGQ finite elements to two-dimensional problems in magnetostatics. It is, essentially, a two-dimensional sequel of step-97. This tutorial discusses two-dimensional curl-curl partial differential equation in four possible configurations. It illustrates the numerical solution to the curl-curl equation in one of the four possible configurations - the two-dimensional vector planar configuration. Experimenting with the gauging parameter, \(\eta^2\), is suggested as a possibility for extension of this tutorial.
The program of this tutorial uses gmsh API. For this reason, the deal.II must be installed with the interface to gmsh. The installation instructions can be found in the README file.
Initially, all problems in electromagnetics are three-dimensional. Reduction of a three-dimensional problem to two dimensions is possible if the initial three-dimensional problem exhibits a translation or rotation symmetry. In both cases the two-dimensional problem domain is a plane in the initial three-dimensional problem domain as shown in the figure below.
There are two types of possible symmetries in the initial three-dimensional problem domain. The magnetic vector potential, \(\vec{A}\), can end up in the two-dimensional problem domain either as an in-plane vector or an out-of-plane vector. Therefore, there are four combinations in total as illustrated by the table below.
| Translation | Rotation | |
| \(\boldsymbol{\vec{A}}\) is an out-of-plane vector | 1) The scalar planar problem | 2) The scalar axisymmetric problem |
| \(\boldsymbol{\vec{A}}\) is an in-plane vector | 3) The vector planar problem | 4) The problem cannot be formulated |
In the scalar planar and scalar axisymmetric problems (numbers 1 and 2 in the table) the curl-curl equation can be replaced with an equivalent div-grad equation. A justification of such replacement can be found in the Replacing curl-curl with div-grad section. The div-grad equation is easier to solve. The magnetic vector potential in these two configurations is gauged by the symmetry in the initial three-dimensional problem domain so no extra effort for gauging is required. Furthermore, there is no need for the current vector potential, \(\vec{T}\), in these two configurations. The compatibility condition and the current vector potential are tools for solving the curl-curl equation, not the div-grad equation.
If the magnetic vector potential is an in-plane vector and the initial three dimensional problem exhibits rotation symmetry (number 4 in the table), the two-dimensional problem cannot be formulated. To see why, let us consider a toroidal inductor. Insert A) in the figure below illustrates the top view of the inductor.
\begin{equation} \vec{\nabla} \cdot \vec{J}_f \ne 0 \end{equation}
in the two-dimensional problem domain. That is, the two-dimensional current density is not a purely solenoidal vector field. On the other hand, the left-hand side of the curl-curl equation,
\[ \vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A}\bigg) = \vec{J}_f, \]
is a curl of a vector field. A curl of a vector field is purely solenoidal. This is true for the two-dimensional vector curl as well, see below. Therefore, the purely solenoidal left-hand side of the curl-curl equation will never match the current density \(\vec{J}_f\) on the right-hand side as \(\vec{J}_f\) is a sum of a solenoidal and a conservative vector fields. It is impossible to solve an equation like that. Or, more precisely, this is not a valid partial differential equation to begin with. In my opinion this unfortunate disposition will persist in any vector axisymmetric problem (number 4 in the table). Consequently, we assume that this type of problems cannot be formulated in terms of the curl-curl partial differential equation.
In this tutorial we will consider the vector planar problem (number 3 in the table). The magnetic vector potential, \(\vec{A}\), in this configuration is an in-plane vector. The initial three-dimensional problem is reduced to two dimensions by exploiting the translation symmetry.
Strictly speaking, the cross product and the curl exist in the three-dimensional space only. To facilitate the analysis we need to introduce two artificial two-dimensional cross products and two artificial curls. The regular cross product in a three-dimensional space takes two vectors as an input. There are two combinations that make sense in two dimensions: (i) both input vectors are in-plane vectors, (ii) one of the input vectors is an in-plane vector while another input vector is an out-of-plane vector. The third combination, i.e., both input vectors are out-of-plane vectors, is not very interesting as a cross product of two collinear vectors equals zero.
The first combination suggests a two-dimensional cross product that takes two in-plane vectors as an input,
\begin{equation} C = \vec{A} \overset{S}{\times} \vec{B} = \begin{vmatrix} 0 & 0 & 1 \\ A_x & A_y & 0 \\ B_x & B_y & 0 \\ \end{vmatrix} = A_x B_y - A_y B_x. \end{equation}
We will call this a two-dimensional scalar cross product as it yields an out-of-plane vector, a scalar. The letter "S" above the cross stands for "scalar".
The second combination suggests a cross product that takes as an input one in-plane vector and one out-of-plane vector,
\begin{equation} \vec{C} = \vec{A} \overset{V}{\times} B = \begin{vmatrix} \hat{i} & \hat{j} & 0 \\ A_x & A_y & 0 \\ 0 & 0 & B \\ \end{vmatrix} = A_y B \hat{i} - A_x B \hat{j}. \end{equation}
We will call this a two-dimensional vector cross product as it yields an in-plane vector. The letter "V" above the cross stands for "vector". There is an alternative definition:
\begin{equation} \vec{C} = B \overset{V}{\times} \vec{A} = \begin{vmatrix} \hat{i} & \hat{j} & 0 \\ 0 & 0 & B \\ A_x & A_y & 0 \\ \end{vmatrix} = - A_y B \hat{i} + A_x B \hat{j}. \end{equation}
From the last two equations we can deduce that the two-dimensional vector cross product is anticommutative, just as the three-dimensional cross product is:
\begin{equation} \vec{A} \overset{V}{\times} B = - B \overset{V}{\times} \vec{A}. \end{equation}
Similarly, we introduce two two-dimensional curls, the scalar curl and the vector curl. First, we assume that a vector field \(\vec{F}\) has the following form in the initial three-dimensional problem domain:
\begin{equation} \vec{F} = F_x(x,y) \hat{i} + F_y(x,y) \hat{j} + 0 \hat{k}. \end{equation}
Then we introduce the scalar curl as
\begin{equation} \vec{\nabla} \overset{S}{\times} \vec{F} = \frac{\partial F_y }{\partial x} - \frac{\partial F_x }{\partial y}. \end{equation}
Next, we assume that a vector field \(\vec{F}\) has the following form in the initial three-dimensional problem domain:
\begin{equation} \vec{F} = 0 \hat{i} + 0 \hat{j} + F(x,y) \hat{k} \end{equation}
and introduce the vector curl as
\begin{equation} \vec{\nabla} \overset{V}{\times} F = \frac{\partial F }{\partial y} \hat{i} - \frac{\partial F }{\partial x}\hat{j}. \end{equation}
The two-dimensional curls retain some properties of the three-dimensional curl. For instance, the curl of a gradient always equals zero,
\begin{equation} \vec{\nabla} \overset{S}{\times} \big( \vec{\nabla} \Phi \big) = \frac{\partial^2 \Phi}{\partial x \partial y} - \frac{\partial^2 \Phi}{\partial x \partial y} = 0. \end{equation}
Similarly, the divergence of a curl always equals zero as well,
\begin{equation} \vec{\nabla} \cdot \big( \vec{\nabla} \overset{V}{\times} F \big) = \frac{\partial^2 F}{\partial x \partial y} - \frac{\partial^2 F}{\partial x \partial y} = 0. \end{equation}
That is to say, the two-dimensional curl yields a purely solenoidal vector field.
Suppose we use the FE_Nedelec finite elements in two dimensions. At a certain point in the program we have decided to ask an object derived from the FEValues class template to compute a curl for us. A good question is: what will we get, the scalar curl or the vector curl? The shape functions of the two-dimensional FE_Nedelec finite elements, \(\vec{\phi}_i\), are in-plane vector fields. Only the scalar curl can operate on in-plane vector fields, see above. Therefore, we will get the scalar curl, \(\vec{\nabla}\overset{S}{\times}\vec{\phi}_i\). The scalar curl yields an out-of-plane vector field, i.e., a scalar field. We must be prepared to get a scalar field when we ask FEValues to compute a curl in two dimensions. The type of curl in deal.II library is defined as Tensor<1, 3> in three dimensions and as a scalar type in two dimensions, see FEValuesViews::Vector::curl_type. The last means that FEValues computes the regular curl \(\big(\vec{\nabla}\times\big)\) in three dimensions and the scalar curl \(\big(\vec{\nabla}\overset{S}\times\big)\) in two dimensions.
If the magnetic vector potential, \(\vec{A}\), is an out-of-plane vector (numbers 1 and 2 in the table above), the curl-curl partial differential equation can be replaced with the div-grad partial differential equation. In this section we will justify this replacement. The information presented in this section is not essential for understanding the rest of the tutorial.
Let us begin with the planar configuration (number 1 in the table above). In this configuration the magnetic vector potential has the following form in the Cartesian coordinate system:
\[ \vec{A} = 0 \hat{i} + 0 \hat{j} + A(x,y) \hat{k}. \]
The free-current density in this configuration has the same form:
\[ \vec{J}_f = 0 \hat{i} + 0 \hat{j} + J_f(x,y) \hat{k}. \]
The last can be verified by a straightforward evaluation of the following expression: \(\vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A} \bigg)\). We also assume that the permeability exhibits the translation symmetry,
\[ \mu = \mu(x,y). \]
The curl of the magnetic vector potential in this configuration can be evaluated as
\begin{equation} \vec{\nabla} \times \vec{A} = \begin{vmatrix} \hat{i} & \hat{j} & \hat{k} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & 0 \\ 0 & 0 & A \\ \end{vmatrix} = \underbrace{\frac{\partial A}{\partial y}}_{\mu F_x} \hat{i} \underbrace{-\frac{\partial A}{\partial x}}_{\mu F_y} \hat{j}. \end{equation}
Next, we introduce the following notation:
\[ \vec{F} = F_x(x,y)\hat{i} + F_y(x,y)\hat{j} = \frac{1}{\mu(x,y)}\bigg( \vec{\nabla} \times \vec{A}(x,y) \bigg). \]
The curl of this vector field can be evaluated as
\begin{equation} \vec{\nabla} \times \vec{F} = \begin{vmatrix} \hat{i} & \hat{j} & \hat{k} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ F_x & F_y & 0 \\ \end{vmatrix} = \bigg[ \frac{\partial F_y}{\partial x} - \frac{\partial F_x}{\partial y} \bigg] \hat{k}. \end{equation}
Next, we express the left-hand side of the three-dimensional curl-curl equation as
\begin{equation} \begin{aligned} \vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A} \bigg) = \vec{\nabla}\times\vec{F} = \bigg[ \frac{\partial F_y}{\partial x} - \frac{\partial F_x}{\partial y} \bigg] \hat{k} = \\ =\bigg[ -\frac{\partial}{\partial x} \frac{1}{\mu} \frac{\partial A}{\partial x} -\frac{\partial}{\partial y} \frac{1}{\mu} \frac{\partial A}{\partial y} \bigg] \hat{k} = -\vec{\nabla} \cdot \bigg( \frac{1}{\mu} \vec{\nabla} A\bigg) \hat{k}. \end{aligned} \end{equation}
Therefore, we can replace the initial three-dimensional curl-curl equation,
\[ \vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A} \bigg) = \vec{J}_f, \]
with
\[ -\vec{\nabla} \cdot \bigg( \frac{1}{\mu} \vec{\nabla} A\bigg) \hat{k} = J_f \hat{k}. \]
Finally, we discard \(\hat{k}\) and arrive at the div-grad equation,
\[ -\vec{\nabla} \cdot \bigg( \frac{1}{\mu} \vec{\nabla} A\bigg) = J_f. \]
In effect, we have replaced the three-dimensional curl-curl equation with the two-dimensional div-grad equation.
Next, let us consider the axisymmetric configuration (number 2 in the table above). We choose to work in the cylindrical coordinate system \((r,\phi,z)\). The reason for this choice is the following. The functionals related to the div-grad equations contain shape functions and their gradients. For this reason, to solve a div-grad equation by finite element method one must be able to evaluate functions and their gradients. The gradient in cylindrical coordinate system under rotation symmetry has the following form
\[ \vec{\nabla}' A(r,z) = \frac{\partial A}{\partial r} + \frac{\partial A}{\partial z}. \]
This is similar to the expression for the gradient in the two-dimensional Cartesian coordinate system,
\[ \vec{\nabla} A(x,y) = \frac{\partial A}{\partial x} + \frac{\partial A}{\partial y}. \]
The function on the \(rz\) plane is evaluated in the exact the same way as on \(xy\) plane. Therefore, for the purpose of solving the div-grad equation by the finite element method we can treat the \(rz\) plane (half-plane, to be precise) as a two-dimensional Cartesian coordinate system with coordinates \(r\) and \(z\). This is an advantage. The spherical coordinate system does not offer this advantage.
In this configuration the magnetic vector potential has the following form in the cylindrical coordinate system:
\[ \vec{A} = 0 \hat{r} + A(r,z) \hat{\phi} + 0 \hat{z}. \]
The free-current density in this configuration has the same form:
\[ \vec{J}_f = 0 \hat{r} + J_f(r,z) \hat{\phi} + 0 \hat{z}. \]
We also assume that the permeability exhibits the rotation symmetry,
\[ \mu = \mu(r,z). \]
The curl of the magnetic vector potential in this configuration can be evaluated as
\[ \vec{\nabla}\times\vec{A} = \underbrace{- \frac{\partial A}{\partial z}}_{\mu F_r} \hat{r} + \underbrace{ \frac{1}{r} \frac{\partial}{\partial r} (r A)}_{\mu F_z} \hat{z}. \]
Next, we introduce the following notation:
\[ \vec{F} = F_r(r,z)\hat{r} + F_z(r,z)\hat{z} = \frac{1}{\mu(r,z)}\bigg( \vec{\nabla} \times \vec{A}(r,z) \bigg). \]
The curl of this vector field in the cylindrical coordinate system is evaluated as
\[ \vec{\nabla}\times\vec{F} = \bigg[ \frac{\partial F_r}{\partial z} -\frac{\partial F_z}{\partial r} \bigg]\hat{\phi} \]
Next, we express the left-hand side of the three-dimensional curl-curl equation equation as
\begin{equation} \begin{aligned} \vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A} \bigg) = \vec{\nabla}\times\vec{F} = \bigg[ \frac{\partial F_r}{\partial z} - \frac{\partial F_z}{\partial r} \bigg] \hat{\phi} = \\ =\bigg[ -\frac{\partial}{\partial z} \frac{1}{\mu} \frac{\partial A}{\partial z} -\frac{\partial}{\partial r} \frac{1}{\mu} \frac{1}{r} \frac{\partial}{\partial r} \big(r A \big) \bigg] \hat{\phi} =\\ =\bigg[ -\frac{\partial}{\partial z} \frac{1}{\mu} \frac{1}{r} \frac{\partial}{\partial z} \big(r A \big) -\frac{\partial}{\partial r} \frac{1}{\mu} \frac{1}{r} \frac{\partial}{\partial r} \big(r A \big) \bigg] \hat{\phi}. \end{aligned} \end{equation}
Then by introducing the following two notations:
\[ \mu' = r \mu \]
and
\[ A' = r A, \]
we can rewrite the left-hand side of the curl-curl equations as
\begin{equation} \begin{aligned} \vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A} \bigg) = \bigg[ -\frac{\partial}{\partial z} \frac{1}{\mu'} \frac{\partial A'}{\partial z} -\frac{\partial}{\partial r} \frac{1}{\mu'} \frac{\partial A'}{\partial r} \bigg] \hat{\phi} = -\vec{\nabla}' \cdot \bigg( \frac{1}{\mu'} \vec{\nabla}' A'\bigg) \hat{\phi}. \end{aligned} \end{equation}
Therefore, we can replace the initial three-dimensional curl-curl equation,
\[ \vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A} \bigg) = \vec{J}_f, \]
with
\[ -\vec{\nabla}' \cdot \bigg( \frac{1}{\mu'} \vec{\nabla}' A'\bigg) \hat{\phi}= J_f\hat{\phi}. \]
Finally, we discard \(\hat{\phi}\) and again arrive at the div-grad equation,
\[ -\vec{\nabla}' \cdot \bigg( \frac{1}{\mu'} \vec{\nabla}' A'\bigg) = J_f. \]
In effect, we have replaced the three-dimensional curl-curl equation in the cylindrical coordinate system with the two-dimensional div-grad equation in the Cartesian coordinate system. Note, that when solving the div-grad equation, the \(r\) and \(z\) coordinates must be treated as coordinates of the two-dimensional Cartesian coordinate system. They must, however, be treated as coordinates of the cylindrical coordinates system when interpreting or postprocessing the results.
The take-away message from this section is that in the two cases of symmetry considered above, where you originally had the curl-curl equation in three dimensions, the formulation of the equation in two dimensions can either use the curl-curl operator or the div-grad operator – the two are mathematically equivalent.
In this tutorial we would like to calculate the magnetic field induced by an infinitely long solenoid equipped with a magnetic core. The figure below illustrates a cross section of the solenoid. The solenoid consists of a magnetic core made of soft magnetic material and current-carrying windings around it. The magnetic core is depicted in green color. The radius of the cross section of the core is \(b_1\). The permeability of the entire space can be described as
\begin{equation} \mu = \left\{ \begin{aligned} & \mu_1 &\text{if } \text{ } & r < b_1 \\ & \mu_0 &\text{if } \text{ } & r > b_1, \end{aligned} \right. \end{equation}
where \(\mu_0\) is the permeability of the free space. The magnetic material is assumed to be homogeneous, linear, lossless, and isotropic.
The blue ring in the figure designates the region filled with current-carrying windings. The current-carrying windings are modeled by a prescribed free-current density. The free-current density can be expressed as
\[ \vec{J}_f = \left\{ \begin{aligned} & 0 &&\text{if }&& r < a_2\\ & K_0 (-y\hat{i} + x\hat{j})&&\text{if }&& a_2 \le r \le b_2\\ & 0 &&\text{if }&& r > b_2\\ \end{aligned} \right. \qquad \qquad \textrm{[Equation 1]} \]
in the Cartesian coordinate system. In the last two equations
\[ r = \sqrt{x^2+y^2}. \]
The magnetic field generated by the solenoid is parallel to the axis of the solenoid. It points perpendicular to the page in the figure above. The magnetic field equals zero outside the solenoid. Note that the word "solenoid" refers to the whole object, the current-carrying windings plus the magnetic core.
The figure below illustrates the corresponding two-dimensional problem domain. The circle \(\Gamma_{I1}\) represents the interface between dissimilar materials. The circle \(\Gamma_{R1}\) represents the outer boundary of the problem domain. We can choose any concentric circle outside the solenoid to be \(\Gamma_{R1}\) as the magnetic field equals zero everywhere outside the solenoid.
The closed-form analytical expression for the magnetic field induced by the solenoid can be easily derived by evaluating Ampere's law,
\[ \int_{\Gamma} \vec{H} = \int_{\Omega} \vec{J}_f, \]
over four rectangular loops \(\Gamma\) drawn in the plane that contains the axis of the solenoid. Example 5.9 in [114] illustrates how to use such Amperian loops. The result reads
\[ \vec{B} = \left\{ \begin{aligned} & 0 &&\text{if }&& r > b_2\\ & \frac{\mu_0 K_0}{2} (b_2^2-r^2)\hat{k}&&\text{if }&& a_2 < r < b_2\\ & \frac{\mu_0 K_0}{2} (b_2^2-a_2^2)\hat{k}&&\text{if }&& b_1 < r < a_2\\ & \frac{\mu_1 K_0}{2} (b_2^2-a_2^2)\hat{k}&&\text{if }&& r < b_1. \end{aligned} \right. \qquad \qquad \textrm{[Equation 2]} \]
The current vector potential of the solenoid has the following form in the three-dimensional Cartesian coordinate system:
\[ \vec{T} = 0\hat{i} + 0\hat{j} + T(x,y)\hat{k}. \]
A straightforward evaluation of the divergence yields
\[ \vec{\nabla} \cdot \vec{T} = 0, \]
which is the Coulomb gauge. That is to say, the current vector potential is gauged by the symmetry of the problem. For this reason, we can predict the closed-form analytical expression for the current vector potential in advance. It is easy to verify that the following expression
\[ \vec{T} = \left\{ \begin{aligned} & \frac{K_0}{2} (b_2^2-a_2^2) \hat{k} &&\text{if }&& r \le a_2 \\ & - \frac{K_0}{2} (x^2+y^2-b_2^2) \hat{k} && \text{if } && a_2 \le r \le b_2 \\ & 0 &&\text{if }&& r \ge b_2, \end{aligned} \right. \qquad \qquad \textrm{[Equation 3]} \]
is gauged by the Coulomb gauge, fits the definition of the current vector potential,
\[ \vec{J}_f = \vec{\nabla} \times \vec{T}, \]
and is continuous on interfaces \(r=a_2\) and \(r=b_2\). Therefore, it is a correct expression. The last requirement, the continuity, can be deduced from the Maxwell's equations. Alternatively, one can have a brief look at the Bossavit's diagrams below and deduce that the \(z\)-component of the current vector potential, \(T\), belongs to the \(H(\text{grad})\) function space. The last implies continuity.
The magnetic vector potential of the solenoid has the following form in the three-dimensional Cartesian coordinate system:
\[ \vec{A} = A_x(x,y)\hat{i} + A_y(x,y)\hat{j} + 0\hat{k}. \]
We cannot infer the divergence of \(\vec{A}\) from this equation. Contrary to the case of the current vector potential, \(\vec{T}\), discussed above, the magnetic vector potential, \(\vec{A}\), is not gauged by the symmetry of the problem. Strictly speaking, it must be gauged differently. In this numerical simulation we will gauge the magnetic vector potential implicitly by adding the \(\eta^2 \vec{A}\) term to the curl-curl equation. The same approach to gauging was taken in step-97. See step-97 for more information about the rationale behind this gauging approach. As soon as we will be using the implicit gauge, the conservative part of the solution, \(\vec{A}\), that will appear at the output of the solver will remain unknown. For this reason, we cannot write down a closed-form analytical expression that would correspond to the numerical solution, \(\vec{A}\).
The description of the operation of the \(\eta^2\) parameter given in step-97 is somewhat simplistic. The Possibilities for extensions section in this tutorial is meant to draw a more granular picture of how the \(\eta^2\) parameter operates. For the experiments suggested in the Possibilities for extensions section it would be handy to have a closed-form analytical expression for \(\vec{A}\) gauged by the Coulomb gauge,
\[ \vec{\nabla}\cdot\vec{A} = 0. \]
Such expression can be derived by integrating the expression for \(\vec{B}\) given above. The result reads
\[ \vec{A} = A_{\phi}(r) \hat{\phi}, \]
in the cylindrical coordinate system and
\[ \vec{A} = A_{\phi}(r) \dfrac{(-y\hat{i}+x\hat{j})}{r} \qquad \qquad \textrm{[Equation 4]} \]
in the Cartesian coordinate system. The scalar field \(A_{\phi}(r)\) is computed separately in each subdomain,
\[ A_{\phi}(r) = \left\{ \begin{aligned} &\dfrac{\mu_0 K_0}{2}\dfrac{b_2}{r}\bigg[\dfrac{1}{4} b_2^3+ \dfrac{1}{b_2}\bigg(a_2(b_2^2-a_2^2)\bigg[\dfrac{1}{2}a_2+\dfrac{1}{2}\dfrac{1}{a_2}b_1^2(\mu_r- 1)\bigg] -\dfrac{1}{2}b_2^2a_2^2 + \dfrac{1}{4}a_2^4\bigg)\bigg] &&\text{if }&& r \ge b_2\\ & \dfrac{\mu_0 K_0}{2}\bigg[\dfrac{1}{2}b_2^2 r -\dfrac{1}{4} r^3+ \dfrac{1}{r}\bigg(a_2(b_2^2-a_2^2)\bigg[\dfrac{1}{2}a_2+\dfrac{1}{2}\dfrac{1}{a_2}b_1^2(\mu_r- 1)\bigg] -\dfrac{1}{2}b_2^2a_2^2 + \dfrac{1}{4}a_2^4\bigg)\bigg] &&\text{if }&& a_2 \le r \le b_2\\ & \dfrac{\mu_0 K_0}{2}(b_2^2-a_2^2)\bigg[\dfrac{1}{2} r + \dfrac{1}{2}\dfrac{1}{r}b_1^2(\mu_r-1) \bigg]&&\text{if }&& b_1 \le r \le a_2\\ & \dfrac{\mu_0 K_0}{2}(b_2^2-a_2^2)\bigg[ \dfrac{1}{2} \mu_r r\bigg] &&\text{if }&& r \le b_1, \end{aligned} \right. \]
where
\[ \mu_r = \frac{\mu_1}{\mu_0}. \]
A straightforward evaluation of divergence yields
\[ \vec{\nabla} \cdot \vec{A} = \vec{\nabla} \cdot \bigg(A_{\phi}(r) \hat{\phi}\bigg) = 0. \]
That is to say, \(\vec{A}\) is gauged by the Coulomb gauge. If the magnetic vector potential has the following form:
\[ \vec{A} = A_{\phi}(r) \hat{\phi}, \]
it relates to the magnetic field as
\[ \hat{k}\cdot\vec{B} = \dfrac{1}{r}\dfrac{\partial} {\partial r}\bigg[r A_{\phi}(r)\bigg]. \]
This expression is just the \(z\) component of the curl expressed in cylindrical coordinate system. It is easy to see that this equality is satisfied by the closed-form analytical expressions for \(\vec{B}\) and \(\vec{A}\) given above. Furthermore, \(A_{\phi}\) is continuous on the interfaces \(r=b_1\), \(r=a_2\), and \(r=b_2\). To summarize, the closed-form analytical expression for \(\vec{A}\) is gauged by the Coulomb gauge, its tangential component is continuous on the interfaces, and it yields correct expressions for the magnetic field if differentiated properly. Therefore, it is a correct expression. Note that the last expression can also be obtained by a straightforward application of the rules of differentiation to the expression for the scalar curl ( \(\vec{\nabla}\overset{S}{\times}\), it is introduced above). That is,
\[ B=\vec{\nabla}\overset{S}{\times}\bigg[ A_{\phi}(r) \dfrac{(-y\hat{i}+x\hat{j})}{r}\bigg] = \dfrac{1}{r}\dfrac{\partial} {\partial r}\bigg[r A_{\phi}(r)\bigg]. \]
It is possible to derive the two-dimensional curl-curl equation, boundary conditions, and interface conditions from the Maxwell's equations. There is, however, a simpler approach. We can use an informal procedure to reduce the three-dimensional equations to two dimensions. Let us reduce the three-dimensional curl-curl equation derived in step-97,
\begin{equation} \vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A}\bigg) + \eta^2 \vec{A} = \vec{\nabla}\times\vec{T}, \end{equation}
as an example. In the two-dimensional vector planar problem the magnetic vector potential, \(\vec{A}\), is an in-plane vector field. Therefore, the curl operator that acts on it must be the scalar curl so we place the letter "S" above the cross product,
\begin{equation} \vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg) + \eta^2 \vec{A} = \vec{\nabla}\times\vec{T}. \end{equation}
The scalar curl yields an out-of-plane vector field. Scaling it by a factor of \(\dfrac{1}{\mu}\) yields an out-of-plane vector field. Therefore, the curl operator that acts on it must be the vector curl. Consequently, we place the letter "V" above the cross product,
\begin{equation} \vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg) + \eta^2 \vec{A} = \vec{\nabla}\times\vec{T}. \end{equation}
We can conclude that the left-hand side in the last equation is an in-plane vector field. Therefore, the right-hand side is an in-plane vector field as well. Consequently, the curl on the right-hand side is the vector curl and the current vector potential has no other choice but to be an out-of-plane vector field,
\begin{equation} \vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg) + \eta^2 \vec{A} = \vec{\nabla}\overset{V}{\times} T. \end{equation}
This is the two-dimensional curl-curl equation we would like to solve. The rationale behind the gauging term, \(\eta^2\vec{A}\), and the current vector potential, \(T\), is discussed in step-97. The same holds in two spatial dimensions.
Other equations of the three-dimensional boundary value problem of step-97 can be reduced to two dimensions by invoking the same simple informal procedure. The result is reads
\begin{equation} \begin{array}{lrcll} \text{ } & \vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg) + \eta^2 \vec{A} = \vec{\nabla}\overset{V}{\times} T & \text{in} & \Omega & \text{(i)},\\ \text{(n)}& \dfrac{1}{\mu}\hat{n}\overset{V}{\times}\bigg(\vec{\nabla}\overset{S}{\times}\vec{A}\bigg)=0 &\text{on} & \Gamma_{R1} & \text{(ii)}, \\ \text{(e)}&\hat{n}\overset{S}{\times}\vec{A}_{+} = \hat{n}\overset{S}{\times}\vec{A}_{-}&\text{on}&\Gamma_{I1}&\text{(iii)},\\ \text{(n)}&\dfrac{1}{\mu}_{+}\hat{n}\overset{V}{\times}\bigg(\vec{\nabla}\overset{S}{\times} \vec{A}_{+} \bigg) - \dfrac{1}{\mu}_{-}\hat{n}\overset{V}{\times} \bigg( \vec{\nabla}\overset{S}{\times} \vec{A}_{-} \bigg) = 0 & \text{on} & \Gamma_{I1}&\text{(iv)}. \end{array} \qquad \qquad \textrm{[Equation 5]} \end{equation}
This time, however, we have only one interface, \(\Gamma_{I1}\), and the homogeneous Neumann boundary condition is used instead of the Robin boundary condition, i.e., \(\gamma = 0\) in equation (ii). The magnetic field outside the solenoid equals zero, \(B=0\). For this reason, the Neumann boundary condition is appropriate,
\begin{equation} \dfrac{1}{\mu}\hat{n}\overset{V}{\times}\bigg(\vec{\nabla}\overset{S}{\times}\vec{A}\bigg)= \dfrac{1}{\mu}\hat{n}\overset{V}{\times} B = 0. \end{equation}
In step-97 the artificial surface \(\Gamma_{R1}\) represents infinity. In this tutorial \(\Gamma_{R1}\) is any circle with a radius \(r>b_2\) concentric with the cross section of the solenoid.
Finally, we note that the magnetic vector potential in current configuration is defined as
\[ B = \vec{\nabla}\overset{S}{\times}\vec{A}. \]
The magnetic field, \(B\), is an out-of-plane vector field. It is, essentially, the \(z\) component of the magnetic field in the initial three-dimensional problem.
As discussed above, the current vector potential in the current configuration is gauged by the symmetry of the problem. Consequently, there is no need in the gauging term, i.e., \(\eta^2 = 0\). Application of the simple informal procedure to the three-dimensional curl-curl equation for \(\vec{T}\) derived in step-97,
\[ \vec{\nabla}\times\bigg(\vec{\nabla}\times\vec{T}\bigg) = \vec{\nabla}\times\vec{J}_f, \]
yields
\[ \vec{\nabla}\overset{S}{\times}\bigg(\vec{\nabla}\overset{V}{\times} T\bigg) = \vec{\nabla}\overset{S}{\times}\vec{J}_f. \]
By substituting the definitions of the scalar and vector curls into the left-hand side of the last equation and rearranging the terms we get
\[ -\vec{\nabla} \cdot \bigg(\vec{\nabla} T\bigg) = \vec{\nabla}\overset{S}{\times}\vec{J}_f. \]
That is, the curl-curl equation in this particular case can be replaced with the div-grad equation. The div-grad equation guarantees the unique solution with a precision to a constant, \(T_0\),
\[ -\vec{\nabla} \cdot \bigg(\vec{\nabla} (T+T_0)\bigg) = -\vec{\nabla} \cdot \bigg(\vec{\nabla} T\bigg). \]
We can choose any constant \(T_0\) we like as the current vector potential has no physical meaning. Only its derivative has the physical meaning of free-current density,
\[ \vec{J}_f = \vec{\nabla}\overset{V} {\times} T. \]
I suggest to chose the constant \(T_0\) such that the current vector potential equals zero on \(\Gamma_{R1}\), i.e., the homogeneous Dirichlet boundary condition. This choice will allow us to neglect the integral \(I_{b3-2}\) in the numerical recipe for magnetic vector potential. We compose the boundary value problem by combining the div-grad equation with the homogeneous Dirichlet boundary condition,
\begin{equation} \begin{array}{rrcll} &-\vec{\nabla}\cdot\bigg( \vec{\nabla} T\bigg) = \vec{\nabla}\overset{S}{\times}\vec{J}_f & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &T = 0 & \text{on} & \Gamma_{R1} & \text{(ii)}. \end{array} \qquad \qquad \textrm{[Equation 6]} \end{equation}
The Dirichlet boundary condition is essential. Minimization of the functional will not enforce it. We will enforce it by constraining the system of linear equations. The current vector potential is computed in the free space. There are no interfaces between dissimilar materials in the free space so no interface conditions are required.
As discussed in step-97, a contemplation of the Bossavit's diagram is the simplest method of assigning a particular type of finite elements to a physical quantity. Let us use the contemplation method in this tutorial as well.
Strictly speaking, the curl exists in the three-dimensional space only. In order to facilitate the analysis we have introduced two artificial curls, see above. That is, in the two-dimensional space the curl dissociates in two curls, the scalar curl and the vector curl. Similarly, the Bossavit's diagram in the two-dimensional space dissociates in two diagrams as illustrated by the figure below.
Let us consider assigning the correct type of finite elements to the magnetic field, \(B\), as an example. First, we contemplate the Bossavit's diagrams and observe that the equation we would like to solve,
\[ \vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg) = \vec{J}_f, \]
belongs to the scalar diagram. It is represented by the red path that links \(\vec{A}\) and \(\vec{J}_f\), see insert A). Second, we observe that the magnetic field \(B\) in the scalar diagram belongs to the \(L_2\) functions space. Third, we look at the table below and conclude that the physical quantities that belong to the \(L_2\) function space are modeled by the FE_DGQ finite elements. Therefore, we need to model the magnetic field, \(B\), by the FE_DGQ finite elements.
This choice makes a lot of sense. Consider the following. In the current problem (number 3 in the table above) the magnetic field is an out-of-plane vector field, i.e., a scalar field. In the table below there are two types of the finite elements that model scalar fields, FE_Q and FE_DGQ. We simply need to choose one of them. The normal component of the magnetic field on interfaces between dissimilar materials is continuous. On the contrary, the tangential component is discontinuous. The first figure in the step-97 tutorial provides an example of how such behavior on interfaces looks like. In the current problem the magnetic field is tangential to the interface \(\Gamma_{I1}\). The tangential component is discontinuous on interfaces. That is, the out-of-plane magnetic field \(B\) we would like to model is discontinuous. Therefore, we need to choose the finite elements that can model discontinuous scalar fields. The FE_DGQ finite elements allow the scalar field they model to be discontinuous on the faces of mesh cells. Then we can get the desired behavior of the magnetic field, \(B\), on the interface \(\Gamma_{I1}\) by constructing the mesh such that the interface is made up of cell faces and by choosing the FE_DGQ finite elements. Note, that the FE_Q finite elements allow no discontinuity in the scalar fields they model.
| Funct. space | Finite elements |
|---|---|
| \(H(\text{grad})\) | FE_Q |
| \(H(\text{curl})\) | FE_Nedelec |
| \(H(\text{div})\) | FE_RaviartThomas |
| \(L_2\) | FE_DGQ |
The same procedure can be used for assigning finite elements to the rest of physical quantities, \(T\), \(\vec{A}\), and \(\vec{J}_f\). The result is summarized in the table below.
| Phys. quant. | Finite elements |
|---|---|
| \(T\) | FE_Q |
| \(\vec{A}\) | FE_Nedelec |
| \(\vec{J}_f\) | FE_RaviartThomas |
| \(B\) | FE_DGQ |
The problem of solving the boundary value problem for the magnetic vector potential can be replaced by the problem of minimizing the following functional:
\begin{equation} \begin{aligned} &F(\vec{A}) = \int_{\Omega}\frac{1}{\mu}\bigg|\vec{\nabla}\overset{S}{\times}\vec{A}\bigg|^2 + \eta^2 \int_{\Omega}\mid\vec{A}\mid^2 -2\int_{\Omega} T \bigg( \vec{\nabla}\overset{S}{\times}\vec{A} \bigg) +2\int_{\Gamma_{R1}} T \bigg(\hat{n}\overset{S}{\times}\vec{A}\bigg). \end{aligned} \end{equation}
The magnetic vector potential, \(\vec{A}\), that minimizes the functional above will satisfy the curl-curl partial differential equation (i), the homogeneous Neumann boundary condition (ii), and the interface condition (iv) of the boundary value problem for \(\vec{A}\), see above. The functional, however, is invariant to the interface condition (iii). This condition is enforced by the choice of finite elements. The FE_Nedelec finite elements guarantee continuity of the tangential component of the magnetic vector potential on the faces of mesh cells. As soon as we build the mesh such that the circle \(\Gamma_{I1}\) is made up of cell faces, the tangential component of \(\vec{A}\) will be continuous, i.e., interface condition (iii) will be satisfied.
We can convert the functional into a system of linear equations by following the procedure discussed in step-97. The resultant system matrix and the right-hand side vector read
\begin{equation} \begin{aligned} & A_{ij} = \underbrace{\int_{\Omega}\frac{1}{\mu} \bigg(\vec{\nabla}\overset{S}{\times}\vec{\phi}_i\bigg) \bigg(\vec{\nabla}\overset{S}{\times}\vec{\phi}_j\bigg)}_{I_{a1}} + \underbrace{\eta^2\int_{\Omega} \vec{\phi}_i \cdot \vec{\phi}_j}_{I_{a3}} \end{aligned} \end{equation}
and
\begin{equation} b_i = \underbrace{\int_{\Omega} T \bigg(\vec{\nabla}\overset{S}{\times}\vec{\phi}_i\bigg) }_{I_{b3-1}} - \underbrace{\int_{\Gamma_{R1}} T \bigg(\hat{n}\overset{S}{\times}\vec{\phi}_i\bigg) }_{I_{b3-2} = 0}. \end{equation}
The vector fields \(\vec{\phi}_i\) in the last two equations are the shape functions of the FE_Nedelec finite elements. The current vector potential, \(T\), is a field function,
\[ T(\vec{r})= \sum_i c_i \phi_i(\vec{r}). \]
The scalar fields \(\phi_i\) in the last equation are the shape functions of the FE_Q finite elements. The degrees of freedom, \(c_i\), are the numerical solution to the boundary value problem for the current vector potential. As soon as we force the current vector potential to zero on the boundary \(\Gamma_{R1}\), i.e., the homogeneous Dirichlet boundary condition (ii) in the boundary value problem for \(T\) above, we can neglect the integral \(I_{b3-2}\) in the recipe for the right-hand side \(b_i\).
The problem of solving the boundary value problem for the current vector potential can be replaced by the problem of minimizing the following functional:
\begin{equation} \begin{aligned} &F(T) = \int_{\Omega}\big|\vec{\nabla} T \big|^2 -2\int_{\Omega} \vec{J}_f \cdot \bigg( \vec{\nabla}\overset{V}{\times}T \bigg) +2\int_{\Gamma_{R1}} \vec{J}_f \cdot \bigg(\hat{n}\overset{V}{\times} T \bigg). \end{aligned} \end{equation}
This functional is invariant to the Dirichlet boundary condition (ii). We will enforce this condition by restricting the degrees of freedom in the system of linear equations. To do so we will invoke the deal.II function VectorTools::interpolate_boundary_values().
The functional above translates into the following numerical recipes for the system matrix and the right-hand side column vector:
\begin{equation} \begin{aligned} & A_{ij} = \underbrace{\int_{\Omega} \bigg(\vec{\nabla} \phi_i\bigg) \cdot \bigg(\vec{\nabla} \phi_j\bigg)}_{I_{a1}} \end{aligned} \end{equation}
and
\begin{equation} b_i = \underbrace{\int_{\Omega}\vec{J}_f\cdot\bigg(\vec{\nabla}\overset{V}{\times} \phi_i\bigg) }_{I_{b3-1}} - \underbrace{\int_{\Gamma_{R1}}\vec{J}_f\cdot\bigg(\hat{n}\overset{V}{\times} \phi_i\bigg) }_{I_{b3-2}=0}. \end{equation}
The scalar fields \(\phi_i\) in the last two equations are the shape functions of the FE_Q finite elements. The free-current density, \(\vec{J}_f\), is defined by the closed-form analytical expression given by the formulation of the problem. As \(\vec{J}_f\) equals zero on \(\Gamma_{R1}\), we can neglect the integral \(I_{b3-2}\) in the recipe for the right-hand side vector \(b_i\).
The free-current density, \(\vec{J}_f\), is computed as
\begin{equation} \vec{J}_f = \vec{\nabla} \overset{V}{\times} T. \qquad \qquad \textrm{[Equation 7]} \end{equation}
This equation seems to be fundamentally different from the partial differential equations. Indeed, to solve a partial differential equation, say
\begin{equation} -\vec{\nabla}\cdot\bigg( \vec{\nabla} T\bigg) = \vec{\nabla}\overset{S}{\times}\vec{J}_f, \qquad \qquad \textrm{[Equation 8]} \end{equation}
one needs to find an unknown function \(T\) that would satisfy the equation. On the contrary, Equation 7 is a direct evaluation of the vector curl. On the higher abstract level, however, there is no difference between Equation 7 and Equation 8. Indeed, we can rewrite Equation 7 as
\[ I\{\vec{J}_f\} = \vec{\nabla} \overset{V}{\times} T, \]
where \(I\{...\}\) is the identity operator. In this form Equation 7 can be perceived as a partial differential equation with the most boring partial differential operator, the identity operator. We can conclude that Equation 7, as many other partial differential equations, can be solved by minimizing a functional. In this particular case, the functional reads
\[ F(\vec{J}_f) = \int_{\Omega} \big| \vec{J}_f \big|^2 - 2\int_{\Omega} \bigg( \vec{\nabla}\overset{V}{\times} T \bigg) \cdot \vec{J}_f. \]
It translates into the following numerical recipes for the system matrix and the right-hand side column vector:
\[ A_{ij} = \underbrace{\int_{\Omega} \vec{\phi}_i \cdot \vec{\phi}_j}_{I_a} \]
and
\[ b_i = \underbrace{\int_{\Omega} \bigg( \vec{\nabla}\overset{V}{\times} T \bigg) \cdot \vec{\phi}_i}_{I_b}. \]
Note that \(A_{ij}\) is the mass matrix. The identity operator always yields the mass matrix. The vector fields \(\vec{\phi}_i\) in the last two equations are the shape functions of the FE_RaviartThomas finite elements. The current vector potential, \(T\), is a field function,
\[ T(\vec{r})= \sum_i c_i \phi_i(\vec{r}). \]
The scalar fields \(\phi_i\) in the last equation are the shape functions of the FE_Q finite elements. The degrees of freedom, \(c_i\), are the numerical solution to the boundary value problem for the current vector potential.
We will call the algorithm that implements this numerical recipe a projector from \(H(\text{grad})\) to \(H(\text{div})\). In the Bossavit's diagram above this projection is represented by the red arrow that links \(T\) in the \(H(\text{grad})\) function space to \(\vec{J}_f\) in the \(H(\text{div})\) function space.
The magnetic field, \(B\), is computed as
\[ B = \vec{\nabla} \overset{S}{\times} \vec{A}. \]
We replace the problem of computing the magnetic field with the problem of minimization of the following functional:
\[ F(B) = \int_{\Omega} \big| B \big|^2 - 2\int_{\Omega} \bigg( \vec{\nabla}\overset{S}{\times} \vec{A} \bigg) B. \]
This functional translates into the following recipes for the system matrix and the right-hand side vector:
\[ A_{ij} = \underbrace{\int_{\Omega} \phi_i \phi_j}_{I_a} \]
and
\[ b_i = \underbrace{\int_{\Omega} \bigg( \vec{\nabla}\overset{S}{\times} \vec{A} \bigg) \phi_i}_{I_b}. \]
The scalar fields \(\phi_i\) in the last two equations are the shape functions of the FE_DGQ finite elements. The magnetic vector potential, \(\vec{A}\), in the last equation is a field function,
\[ \vec{A}(\vec{r}) = \sum_i c_i \vec{\phi}_i(\vec{r}), \]
where vector fields \(\vec{\phi}_i\) are the shape functions of the FE_Nedelec finite elements. The degrees of freedom \(c_i\) are the numerical solution of the boundary value problem for the magnetic vector potential, see above.
We will call the algorithm that implements this numerical recipe a projector from \(H(\text{curl})\) to \(L_2\). In the Bossavit's diagram above this projection is represented by the red arrow that links \(\vec{A}\) in the \(H(\text{curl})\) function space to \(B\) in the \(L_2\) function space.
Coming back to the concrete problem of solving for the fields that describe the solenoid discussed above, let us now talk about the concrete set-up of the problem. The coarse mesh is constructed with a help of gmsh. The circle.geo file contains the description of the mesh. The mesh file, circle.msh, can be generated by executing the following command.
The step-98 program loads circle.msh and refines the mesh globally r times, where r is the mesh refinement parameter.
The mesh is a discretized version of the problem domain shown in the fourth figure on this page. The mesh is constructed such that all circles in the problem domain are delineated by cell faces. That is, no circle runs through a mesh cell. The following figure illustrates the mesh with r=2.
The main elements of the mesh are listed below.
In this tutorial program all IDs (boundary, material, and manifold) are assigned in the circle.geo file. The program just loads the mesh together with the IDs. Note also that we will use the second version of the function GridIn::read_msh as we would like to pass all IDs to deal.II via physical names, not via physical tags. For example, the circle.geo file contains the following lines.
Here "BoundaryID: 1" is the physical name and 2 next to it is the physical tag. The faces that have elementary tags listed in the curly brackets will be grouped into a singe boundary with boundary ID=1. The physical tag 2 will be discarded.
The table below lists the manifold IDs assigned to various constituents of the mesh and the corresponding manifolds.
| Manifold ID | Manifold |
|---|---|
| 1 | SphericalManifold |
| 2 | FlatManifold |
| 3 | TransfiniteInterpolationManifold |
The next table lists the material IDs and the corresponding materials.
| Material ID | Material |
|---|---|
| 1 | Free space, \(\mu=\mu_0\), \(\vec{J}_f = 0\). |
| 2 | Magnetic material, \(\mu=\mu_1\), \(\vec{J}_f = 0\). |
| 3 | Current-carrying windings of the coil, \(\mu=\mu_0\), \(\vec{J}_f \ne 0\). |
The boundary ID=-1 is reserved for internal boundaries (interfaces between dissimilar materials and interfaced between regions of dissimilar geometries), see the documentation of the second version of the GridIn::read_msh function.
The following table summarizes the boundary IDs.
| Boundary ID | Boundary |
|---|---|
| 1 | The outer boundary of the mesh |
| -1 | All circular interfaces between dissimilar materials and boundary of the square region in the middle of the mesh |
The figure below illustrates the geometry of the mesh together with the assigned IDs.
The assignment of the boundary and material ID is quite straightforward. The assignment of the manifold IDs deserves an explanation. All circular boundaries and interfaces are marked with the manifold ID=1 what corresponds to the spherical manifold. Only the first global refinement in deal.II will apply this manifold, see the manifold topic. As soon as we will refine globally more than one time, we attach the spherical manifold (ID=1) to the surfaces between the circular boundaries and interfaces. This will allow the faces of the cells created by the refinement process to inherit the manifold ID=1. This way cells with circular faces will be created everywhere in between the circular boundaries. The square boundary in the middle of the mesh and the area inside it receive the manifold ID=2 what corresponds to the flat manifold. The flat manifold is appropriate for rectangular cells, see, the documentation of the function GridGenerator::plate_with_a_hole, instance. The area in between the square and the innermost circle is attached to the transfinite interpolation manifold (ID=3). It allows us to accurately fuse the rectangular cells with the cell that have circular faces, see the documentation of the TransfiniteInterpolationManifold class template.
The program uses four meshes of different degrees of refinement. Note that only globally refined meshes are used. There are no non-conforming cells, hanging nodes, and hanging node constrains in this tutorial. We will use constrains only for the purpose of enforcing the Dirichlet boundary condition and for distributing the components of the local (cell specific) system matrices and right-hand sides to the global matrix and right-hand side by calling the function AffineConstraints::distribute_local_to_global.
The program runs in a loop. In each iteration of the loop a different refinement parameter, r, is assumed. Each iteration consists of four stages:
\[ -\vec{\nabla}\cdot\bigg( \vec{\nabla} T\bigg) = \vec{\nabla}\overset{S}{\times}\vec{J}_f. \]
The source on the right-hand side of this equation, \(\vec{J}_f\), is given as a closed-form analytical expression by Equation 1. The quality of the computed current vector potential, \(T\), is assessed by computing the \(L^2\) error norm. The \(L^2\) error norm is computed by comparing the numerical result with the closed-form analytical expression given by Equation 3. The \(L^2\) error norm is saved in a convergence tabletable_T.\[ \vec{J}_f = \vec{\nabla}\overset{V} {\times} T. \]
This is a stand alone equation as there are no boundary and interface conditions to be observed. For this reason, this equation is not associated with a boundary value problem. The quality of the computed free-current density, \(\vec{J}_f\), is assessed by computing the \(L^2\) error norm. The \(L^2\) error norm is computed by comparing the numerical result with the closed-form analytical expression given by Equation 1. The \(L^2\) error norm is saved in a convergence tabletable_Jf. Strictly speaking, the conversion of \(T\) back to \(\vec{J}_f\) is not necessary for achieving the final goal (which is computing the magnetic field). This conversion is done as an extra check of the current vector potential, \(T\), computed at the preceding stage.\[ \vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg) + \eta^2 \vec{A} = \vec{\nabla}\overset{V}{\times} T. \]
The source on the right-hand side of this equation, \(T\), is computed numerically at the first stage above. Recall that we apply an implicit gauge. The last means that we cannot predict in advance the conservative part of the numerical solution, \(\vec{A}\), which will appear at the output of the solver. Consequently, we cannot write down a closed-form analytical expression that will predict the numerical solution. This, in turn, implies that we cannot assess the quality of the numerical solution directly by computing the \(L^2\) error norm. Instead, we will assess the quality indirectly by observing the \(L^2\) error norm computed for the magnetic field at the next stage. Despite the fact that computing \(L^2\) error norms at this stage is meaningless due to the implicit gauge, we will compute them anyway. The reason for this is the following. The Possibilities for extensions section offers a numerical experiment the goal of which is to find a value of \(\eta^2\) which corresponds to Coulomb gauge ( \(\vec{\nabla} \cdot \vec{A} = 0\)). The conservative part of the numerical solution in this case is quite predictable - it equals zero. Consequently, we can write down the closed-form analytical expression that predicts the numerical result in this particular case. Equation 4 gives such an expression. So, at this stage we compute the \(L^2\) error norms for \(\vec{A}\) by comparing the numerical result with Equation 4 and save them intable_A. These error norms can be useful when conducting experiments with the Coulomb gauge as suggested in the Possibilities for extensions section. Keep in mind that in the default configuration, \(\eta^2=0\), table_A contains no useful information.\[ B = \vec{\nabla} \overset{S}{\times} \vec{A}. \]
This is a stand alone equation as there are no boundary and interface conditions to be observed. For this reason, this equation is not associated with a boundary value problem. The quality of the computed magnetic field, \(B\), is assessed by computing the \(L^2\) error norm. The \(L^2\) error norm is computed by comparing the numerical result with the closed-form analytical expression given by Equation 2. The \(L^2\) error norm is saved in a convergence tabletable_B.Each stage is implemented by the code contained in a dedicated namespace. The correspondence between the stages and the namespaces is given in the table below.
| Stage | Namespace |
|---|---|
| 1 | SolverT |
| 2 | ProjectorHgradToHdiv |
| 3 | SolverA |
| 4 | ProjectorHcurlToL2 |
Each of the four namespaces contains a solver which solves for a potential ( \(T\) or \(\vec{A}\)) or converts a potential into a measurable field ( \(\vec{J}_f\) or \(B\)). All four solvers have a lot in common. Each solver setups the simulation, solves the system of linear equations, saves the results, etc. The code common for all solvers is implemented in the BaseSolver class template. The BaseSolver class along with the auxiliary structure UpdateFlagsCollection is implemented in the namespace BaseClasses.
The namespace ExactSolutions contains exact closed-form analytical expressions for \(\vec{J}_f\) (Equation 1), \(B\) (Equation 2), \(T\) (Equation 3), and \(\vec{A}\) (Equation 4).
The following is the control panel of the program. The scaling of the program can be changed by setting mu_0 = 1.0. Then the computed magnetic field, \(B\), will have to be multiplied by a factor of \(\mu_0 = 1.25664 \cdot 10^{-6}\). The free-current density, \(\vec{J}_f\), and the current vector potential, \(T\), do not depend on scaling. The parameter fe_degree encodes the degree of the FE_Nedelec, FE_RaviartThomas, and FE_DGQ finite elements. The degree of the FE_Q finite elements is computed as fe_degree + 1. The reason for this is that the lowermost degree of the FE_Q finite elements equals 1 while the lowermost degree of the FE_Nedelec, FE_RaviartThomas, and FE_DGQ finite elements equals 0 by convention. The boundary, material, and manifold IDs are set in the circle.geo file. The IDs listed below must match the corresponding IDs set in the circle.geo file. If project_exact_solution = true, the program projects the exact solutions for \(T\), \(\vec{J}_f\), \(\vec{A}\), and \(B\) onto the corresponding function spaces and saves the results into the corresponding .vtu files next to the numerical solutions. By default, this feature is switched off. It is used only in the Possibilities for extensions section.
The following class describes a convergence table. The convergence tables are saved on disk in TeX format.
The following namespace contains closed-form analytical expressions for \(T\), \(\vec{J}_f\), \(\vec{A}\), \(B\), mentioned in the introduction to this tutorial.
The following function describes the free-current density, \(\vec{J}_f\), inside the current region. The current density in this tutorial is implemented by ExactSolutions::FreeCurrentDensity class and by SolverT::Solver::free_current_density member function . There is a subtle difference in how these two implementation compute the free-current density. Both classes, however, utilize the same expression for the free-current density. This function describes the expression.
The following class implements the closed-form analytical expression for the free-current density, \(\vec{J}_f\), in the entire domain. The free-current density is computed purely on the basis of the spatial coordinates of the field point. In go coordinates, out comes the free-current density. The information on the material ID of the mesh cells and any other information on the mesh is ignored. This function is used for computing \(L_2\) error norms and for computing the projected exact solution. The \(\vec{J}_f\) on the right-hand side of the div-grad equation is implemented by the member function SolverT::Solver::free_current_density
The following class implements the closed-form analytical expression for the magnetic vector potential, \(\vec{A}\).
The following class implements the closed-form analytical expression for the current vector potential, \(T\).
The following class implements the closed-form analytical expression for the magnetic field, \(B\).
As discussed above, the following namespace aggregates the code common to all four solvers used in the tutorial. All four solvers are derived from the BaseSolver class.
The computation of the integrals of the functionals is delegated to the derived classes. The function that computes the integrals, BaseSolver::system_matrix_local, is virtual and must be overridden by the derived classes. However, the objects of the types FEValues are initialized at the level of the BaseSolver class together with AssemblyScratchData. For the initialization to work properly, the BaseSolver class must know which cell data to compute for a particular implementation of the solver down the hierarchy. This information is communicated to the BaseSolver by passing an argument of the type UpdateFlagsCollection to the constructor. In all solvers, with exception of SolverT, we have two types of finite elements. One type of finite elements models the solution to the partial differential equation. Another models the physical quantity on the right-hand side of the partial differential equation. The solution_update_flags data member below contains the flags for updating the values of the finite elements that model the solution. The rhs_update_flags data member contains the flags for updating the values of the finite elements that model the physical quantity on the right-hand side of the partial differential equation.
Each iteration of the program consists of four stages. Each stage utilizes one solver. All four solvers used in this tutorial are derived from the following class. The solver used in the first stage loads the mesh and refines it if necessary. The solvers in the other three stages reuse the mesh prepared at the first stage. Furthermore, the solver in the first stage expects a closed-form analytical expression on the right-hand side of the partial differential equation. Each solver in the other three stages expects a potential computed at one of the preceding stages, i.e., a field in a form of linear superposition of the shape functions. At the top of the following class, we declare two constructors. The first constructor must be used for constructing solver for the first stage. The second constructor must be used for constructing the solvers for the second, third, and fourth stages. The second constructor has three extra arguments, triangulation_rhs, dof_handler_rhs, and solution_rhs for accommodating the mesh and the potential computed at one of the preceding stages.
Following the constructors are the eight functions that implement various steps typical for every solver. The function run aggregates these steps. This arrangement of functions is quite standard in deal.II, see step-3, for instance. Normally, the setup function begins by distributing the dofs. The problem is: the type of the finite elements is not known in the BaseSolver class. It is specified in the derived classes. Therefore, it could be reasonable to make the setup function virtual as well. Instead, we move the dof distribution code to the end of the make_mesh function which is, in fact, virtual and must be overridden in the derived classes anyway.
After that, we declare six get functions that provide access to protected data. These functions are called from outside the solver, see MagneticProblem::run function. The data provided by these get-functions is used to fill the convergence tables and to pass the references to the triangulation, the dof handler, and the dofs between solvers.
We begin the protected section of the BaseSolver class by declaring three data members that store the input from one of the preceding solvers. If there is no preceding solver and these data is not provided, i.e., first of the two constructors above has been used, these three data members point to triangulation, dof_handler, and solution, of the current solver.
Following are the three data members that store the triangulation, dof handler, and dofs vector of the current solver. In the case of the first-stage solver, the triangulation data member stores the loaded and refined mesh. This data member is not used in the case of the solvers of the second, third, and fourth stages. The dof_handler and solution represent the result of the solver, i.e., the computed potential or field.
Next, we declare four data members that describe the system of linear equations to be solved by the linear solver. The components of the system matrix and that of the right-hand side are computed by the assemble function. The affine constraints are used to apply the Dirichlet boundary conditions and to distribute the local (cell specific) system matrix and right-hand side to system_matrix and system_rhs. The are no hanging nodes and hanging node constraints in this program. The last data member in this block describes the dynamic sparsity pattern. The tutorial step-2 discusses the rationale behind the dynamic sparsity pattern.
The next block contains two data members related to the exact solution. The data member exact_solution points to the closed-form analytical solution the solver attempts to compute. If Settings::project_exact_solution=true, the exact solution is projected onto a proper function space and the data member projected_exact_solution is populated by the dofs of the projected exact solution. The corresponding dof handler is dof_handler. Together dof_handler and projected_exact_solution constitute the field function which describes the exact solution. It is saved into the .vtu file next to the solution and the \(L_2\) error norm.
The next four data members simply store the data supplied as arguments to the constructor. The stage data member stores the number of the current stage. The mapping_degree data member contains the degree of mapping from the reference cell to a mesh cell and back. The update_flags_collection contains the information on which finite element values must be computed for each cell. The names of the output files are derived by appending strings to file_name.
The next block contains two data members computed by the function BaseSolver::compute_error_norms. The L2_per_cell data member contains one value of the \(L^2\) error norm per mesh cell. It is saved into the .vtu file next to the solution. The L2_norm data member contains one value of the \(L^2\) error norm per mesh. It is reported in the convergence table.
The program utilizes the WorkStream technology. The step-9 tutorial does a much better job of explaining the workings of WorkStream. Reading the WorkStream paper is recommended. In very simple terms, the workings of the WorkStream can be envisioned as the following. Let us assume we have a task of computing components of the system matrix, \(A_{ij}\), and the right-hand side vector, \(b_i\). Simply put, we need to fill in the matrix system_matrix and vector system_rhs. This is a big task as the number of dofs is large. The idea is to split the big task on a number of small tasks and feed them to multiple threads to speed up the calculation process. Each small task consists of computing the contributions of a single mesh cell to system_matrix and system_rhs. These contributions are stored temporary in cell_matrix and cell_rhs for each cell. These contributions are then copied to system_matrix and system_rhs. WorkStream creates and schedules the small tasks and takes care of copying cell_matrix and cell_rhs to system_matrix and system_rhs. The rest of the declarations in the protected section of the BaseSolver class help to communicate to WorkStream information which is necessary for its operation.
First, the CellIteratorPair type is declared in the block of code below. The solvers in the second, third, and fourth stages use two dof handlers, dof_handler and dof_handler_rhs. The WorkStream needs to walk through the two dof handlers synchronously. For this purpose we pair two active cell iterators (one from dof_handler, another from dof_handler_rhs). For that we need the CellIteratorPair type. The solver at the first stage (a solver constructed by invoking the first constructor above) uses only one dof handler. In this case the constructor makes dof_handler_rhs to reference dof_handler. In effect, both iterators of the tuple will iterate the same dof handler. In the case of the first-stage solver we will use only the first iterator.
Next, we declare the AssemblyScratchData type. An object of this type contains the relevant information on the current mesh cell which is used as an input for computing the components of the system matrix and the right-hand side. WorkStream creates an object of this type and passes it to the function system_matrix_local which, in turn, computes components of cell_matrix and cell_rhs.
Next, the type AssemblyCopyData is declared. An objects of this type contains the cell specific contributions to the system matrix and the right-hand side. The WorkStream creates object of this type and passes it to function system_matrix_local along with an object of the type AssemblyScratchData. The function system_matrix_local, in turn, takes input data form the object of the type AssemblyScratchData, computes the relevant integrals and places the result into the object of the type AssemblyCopyData. The WorkStream copies the content of the AssemblyCopyData object, cell_matrix and cell_rhs, into system_matrix and system_rhs. The data member AssemblyCopyData::local_dof_indices contains the indices of global components to which cell-specific data must be copied.
WorkStream calls the function system_matrix_local to compute the cell-specific components. Likewise, WorkStream uses calls to copy_local_to_global to copy the cell-specific data into the system matrix and right-hand side. The functions system_matrix_local and copy_local_to_global are declared last.
The following are the implementations of the two constructors of the BaseSolver class.
The following function applies the Dirichlet boundary condition, sets up a sparsity pattern, and initializes the vectors and matrices. It is common for the setup function to distribute the dofs. The type of the finite elements, however, is not known at the level of BaseSolver class. The type of the finite elements is chosen in the derived classes. For this reason, the task of distributing the dofs is shifted to the end of the make_mesh function which is a virtual function.
The Dirichlet boundary condition must be enforced only in the first stage where the div-grad equation is solved for the current vector potential, T. For this reason we have if (stage == 1) filter in the beginning of the function. The boundary value problem for the curl-curl equation utilizes the Neumann boundary condition. It is a natural boundary condition. It is enforced by minimization of the functional. The two projectors, \(T \rightarrow \vec{J}_f\) and \(\vec{A} \rightarrow B\), use no boundary conditions.
Formally, the following function assembles the system of linear equations. In reality, however, it just spells all the magic words to get the WorkStream going. The interesting part, i.e., computing the components of the system matrix and the right-hand side, happens in the function system_matrix_local which is a virtual function. That is to say, the recipe for the functional is implemented in the derived class by overriding the virtual function system_matrix_local.
The following are the implementation of the constructors of the AssemblyScratchData. The first constructor initializes the scratch data from the input parameters. The second - from another object of the same type, i.e., a copy constructor.
The following function copies the components of a cell matrix and a cell right-hand side into the system matrix, \(A_{ij}\), and the system right-hand side, \(b_i\).
The following function solves the system of linear equations. The stopping condition for the iteration algorithm is \(\|\boldsymbol{b}-\boldsymbol{A}\boldsymbol{c}\|<10^{-8}\|\boldsymbol{b}\|\). The maximum number of iteration steps is set to system_rhs.size() as the conjugate gradient algorithm is supposed to find the solution in at most \(m\) steps for an \(m \times m\) system matrix. This function also distributes constraints. The constraints are only used to enforce the Dirichlet boundary condition.
The following two functions compute the error norms and project the exact solution.
The following function saves the solution and the \(L_2\) error norm into a .vtu file. If Settings::project_exact_solution = true, the projected exact solution is saved as well.
The following function clears the memory for the next solver.
The following functions calls all the constituent functions in the right order.
The following is the straightforward implementation of the get functions.
The following namespace contains the code related to the computation of the current vector potential, \(T\).
We derive the solver from the BaseSolver class. What is left to do is to initialize the BaseSolver, override the two virtual functions (make_mesh and system_matrix_local), and implement the function free_current_density. The function free_current_density implements the closed-form analytical expression for \(\vec{J}_f\) on the right-hand side of the div-grad equation.
Following is the implementation of the constructor. We use the first constructor of the BaseSolver class as the solver is used at the first stage. By looking at the expressions for \(A_{ij}\) and \(b_i\) we can conclude that to compute them we need gradients of the shape functions, quadrature points, and the quadrature weights multiplied by the Jacobian determinant(JxW). The quadrature points are needed to sample the closed-form analytical expression for \(\vec{J}_f\). Accordingly, we use the update flags update_quadrature_points, update_gradients, and update_JxW_values for the FE_Q finite elements. This time we do not use the finite elements that model the potential on the right-hand side of the equation, so we set update_default for the right-hand side finite elements.
The following function loads the mesh, creates manifolds, bounds the manifolds to the manifold IDs, refines the mesh, and distributes the dofs. Note that we are allowed to create the manifolds locally in the function as the triangulation object keeps copies of the manifolds, see step-65. Also recall that we have shifted the task of distributing the dofs from the setup function to the make_mesh function to evade the necessity to make the setup function virtual. This allows us to keep one setup function in the BaseSolver class that serves the needs of all derived classes.
The following function assembles a fraction of the system matrix and the system right-hand side related to a single cell. These fractions are copy_data.cell_matrix and copy_data.cell_rhs. They are copied to system_matrix and system_rhs by WorkStream.
The following function implements the closed form analytical expression for \(\vec{J}_f\) on the right-hand side of the div-grad equation.
The following namespace contains all the code related to the computation of the free-current density, \(\vec{J}_f\).
We derive the solver from the BaseSolver class. What is left to do is to initialize the BaseSolver and override the two virtual functions (make_mesh and system_matrix_local).
Following is the implementation of the constructor. We use the second constructor of the BaseSolver class as the solver is used at the second stage. By looking at the expressions for \(A_{ij}\) and \(b_i\) we can conclude that to compute them we need values of the shape functions and the quadrature weights multiplied by the Jacobian determinant(JxW) from the FE_RaviartThomas finite elements. Accordingly, we use the update flags update_values, and update_JxW_values for the FE_RaviartThomas finite elements. This time there is a numerically computed potential, \(T\), on the right-hand side of the equation. It is modeled by the FE_Q finite elements. To compute the right-hand side, we need gradients of the shape functions. Accordingly, we use the update flag update_gradients for the FE_Q finite elements.
At the second stage we do not load the mesh. We reuse the mesh loaded at the first stage. Consequently, we just need to distribute the dofs.
The following function assembles a fraction of the system matrix and the system right-hand side related to a single cell. These fractions are copy_data.cell_matrix and copy_data.cell_rhs. They are copied to system_matrix and system_rhs by WorkStream.
The following namespace contains all the code related to the computation of the magnetic vector potential, \(\vec{A}\).
We derive the solver from the BaseSolver class. What is left to do is to initialize the BaseSolver, override the two virtual functions (make_mesh and system_matrix_local), and implement the function permeability.
Following is the implementation of the constructor. We use the second constructor of the BaseSolver class as the solver is used at the third stage. By looking at the expressions for \(A_{ij}\) and \(b_i\) we can conclude that to compute them we need values of the shape functions, their gradients, and the quadrature weights multiplied by the Jacobian determinant(JxW) from the FE_Nedelec finite elements. Accordingly, we use the update flags update_values, update_gradients, and update_JxW_values for the FE_Nedelec finite elements. This time there is a numerically computed potential, \(T\), on the right-hand side of the equation. It is modeled by the FE_Q finite elements. To compute it, we need values of the shape functions. Accordingly, we use the update flag update_values for the FE_Q finite elements.
At the third stage we do not load the mesh. We reuse the mesh loaded at the first stage. Consequently, we just need to distribute the dofs.
The following function assembles a fraction of the system matrix and the system right-hand side related to a single cell. These fractions are copy_data.cell_matrix and copy_data.cell_rhs. They are copied to system_matrix and system_rhs by WorkStream.
The following function implements the equation for permeability.
The following namespace contains all the code related to the computation of the magnetic field, \(B\).
We derive the solver from the BaseSolver class. What is left to do is to initialize the BaseSolver and override the two virtual functions (make_mesh and system_matrix_local).
Following is the implementation of the constructor. We use the second constructor of the BaseSolver class as the solver is used at the fourth stage. By looking at the expressions for \(A_{ij}\) and \(b_i\) we can conclude that to compute them we need values of the shape functions and the quadrature weights multiplied by the Jacobian determinant(JxW) from the FE_DGQ finite elements. Accordingly, we use the update flags update_values, and update_JxW_values for the FE_DGQ finite elements. This time there is a numerically computed potential, \(\vec{A}\), on the right-hand side of the equation. It is modeled by the FE_Nedelec finite elements. To compute the right-hand side, we need gradients of the shape functions. Accordingly, we use the update flag update_gradients for the FE_Nedelec finite elements.
At the fourth stage we do not load the mesh. We reuse the mesh loaded at the first stage. Consequently, we just need to distribute the dofs.
The following function assembles a fraction of the system matrix and the system right-hand side related to a single cell. These fractions are copy_data.cell_matrix and copy_data.cell_rhs. They are copied to system_matrix and system_rhs by WorkStream.
The MagneticProblem class hosts the main loop inside the run function. The implementation of the loop is straightforward - we create and run the four solvers one-by-one and copy the relevant data into convergence tables.
Stage 1. Computing \(T\).
Stage 2. Computing \(\vec{J}_f\).
Stage 3. Computing \(\vec{A}\).
Stage 4. Computing \(B\).
End stage 4.
gmsh has global state, so set that up in the normal way:
The program generates the following output in the command line interface by default.
The program assumes the finite elements of the lowermost degree, i.e., \(p' = 1\) for the FE_Q finite elements and \(p=0\) for other finite elements. To change the degree of the finite elements, say \(p' = 3\) and \(p = 2\), one needs to change the setting Settings::fe_degree = 2 and rebuild the program. The degree of the FE_Q finite elements will be computed automatically as \(p'= p + 1\).
The program dumps a number of files in the current directory. In the default configuration these files are:
.vtu files. They contain the computed vector fields. Recall that the spherical manifold and transfinite interpolation manifold are attached to many cell faces. Consequently, these cell faces are curved. Furthermore, the shape functions are mapped from the reference cell to the real mesh cells by the second-order mapping to accommodate the cells with curved faces. For these reasons, one needs to use a visualization software that can deal with curved faces and the higher-order mapping. A fresh version of ParaView is recommended. Visit did not have this feature at the time this tutorial was written (in early 2026). The Notes on visualizing high order output provide more information on this topic..tex files. These files contain the convergence tables.The following provides examples of the convergence tables simulated with the default settings for three different degrees of the finite elements, \(p = 0, 1, 2\) (recall that \(p' = p + 1\)).
| p' | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 1 | 1 | 144 | 153 | 7.47e-04 | - |
| 1 | 2 | 576 | 593 | 1.86e-04 | 2.01 |
| 1 | 3 | 2304 | 2337 | 4.64e-05 | 2.00 |
| 1 | 4 | 9216 | 9281 | 1.16e-05 | 2.00 |
| 2 | 1 | 144 | 593 | 1.38e-05 | - |
| 2 | 2 | 576 | 2337 | 8.72e-07 | 3.98 |
| 2 | 3 | 2304 | 9281 | 5.48e-08 | 3.99 |
| 2 | 4 | 9216 | 36993 | 3.45e-09 | 3.99 |
| 3 | 1 | 144 | 1321 | 1.38e-05 | - |
| 3 | 2 | 576 | 5233 | 8.72e-07 | 3.98 |
| 3 | 3 | 2304 | 20833 | 5.48e-08 | 3.99 |
| 3 | 4 | 9216 | 83137 | 3.51e-09 | 3.97 |
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 0 | 1 | 144 | 296 | 2.50e-02 | - |
| 0 | 2 | 576 | 1168 | 1.25e-02 | 1.00 |
| 0 | 3 | 2304 | 4640 | 6.27e-03 | 1.00 |
| 0 | 4 | 9216 | 18496 | 3.13e-03 | 1.00 |
| 1 | 1 | 144 | 1168 | 3.04e-04 | - |
| 1 | 2 | 576 | 4640 | 3.80e-05 | 3.00 |
| 1 | 3 | 2304 | 18496 | 4.74e-06 | 3.00 |
| 1 | 4 | 9216 | 73856 | 5.91e-07 | 3.00 |
| 2 | 1 | 144 | 2616 | 3.04e-04 | - |
| 2 | 2 | 576 | 10416 | 3.79e-05 | 3.00 |
| 2 | 3 | 2304 | 41568 | 4.74e-06 | 3.00 |
| 2 | 4 | 9216 | 166080 | 5.91e-07 | 3.00 |
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 0 | 1 | 144 | 144 | 2.01e-08 | - |
| 0 | 2 | 576 | 576 | 9.72e-09 | 1.05 |
| 0 | 3 | 2304 | 2304 | 4.81e-09 | 1.02 |
| 0 | 4 | 9216 | 9216 | 2.40e-09 | 1.00 |
| 1 | 1 | 144 | 576 | 4.15e-10 | - |
| 1 | 2 | 576 | 2304 | 1.02e-10 | 2.02 |
| 1 | 3 | 2304 | 9216 | 2.54e-11 | 2.01 |
| 1 | 4 | 9216 | 36864 | 6.36e-12 | 2.00 |
| 2 | 1 | 144 | 1296 | 2.32e-11 | - |
| 2 | 2 | 576 | 5184 | 1.46e-12 | 3.99 |
| 2 | 3 | 2304 | 20736 | 9.15e-14 | 3.99 |
| 2 | 4 | 9216 | 82944 | 5.94e-15 | 3.95 |
The following notations were used in the headers of the tables:
Let us contemplate these convergence tables for a brief moment. The first table illustrates convergence of the numerically computed current vector potential, \(T\). In this particular case we can expect the order of the convergence rate of \(\alpha_{L^2} \le p' + 1\) (See also video lecture 3.95.) If the finite elements of the lowermost degree are used, \(p'=1\), the order of convergence rate is at the upper boundary of the expected values, \(\alpha_{L^2} \approx 2.0\). If the finite elements of the second degree are used, \(p'=2\), the order of the convergence rate is higher than expected, \(\alpha_{L^2} \approx 4.0\). This means that the \(L^2\) error norm converges to zero at a rate higher than theoretically possible. Most likely this is due to the fact that the current vector potential has a particularly simple form, Equation 3. It is either constant or changes as the second-order monomial,
\[ T \sim r^2. \]
The second order polynomial approximates this behavior exactly. That is to say, in this particular case we have a lucky situation in which the shape functions can approximate the numerical solution exactly within each mesh cell. This, most likely, explains the extra rapid convergence of the \(L^2\) norm. Note also that the error norms, \(\|e\|_{L^2}\), at \(p'=2\) and \(p'=3\) are the same for the same values of the mesh refinement parameter, \(r\). This is, most likely, due to the fact that the second-order shape functions, \(p'=2\), approximate the solution exactly within each mesh cell and the third-order monomials of the shape function at \(p'=3\) have absolutely nothing to contribute to the quality of approximation.
The second table illustrates convergence of the numerically computed free-current density, \(\vec{J}_f\). The free-current density has been computed as a derivative of current vector potential,
\[ \vec{J}_f = \vec{\nabla}\overset{V} {\times} T. \]
The derivative reduces the order of the convergence rate by one. Most likely in this particular case the error made in computing \(\vec{J}_f\) is defined by the error made in computing \(T\). For this reason, the order of convergence in the second table equals the order of convergence in the first table minus one.
Due to the implicit gauge we cannot observe the convergence of the magnetic vector potential, \(\vec{A}\). Instead, we can observe the convergence on the magnetic field, \(B\), given in the third table. It follows the same pattern: The rate of convergence at the lower degrees of the finite elements is at the best expected value, \(\alpha_{L^2} = p + 1\); the rate of convergence at the higher degrees of the finite elements is better than expected. The extra rapid convergence at the higher degrees of the finite elements is explained by the relatively simple form of the field being approximated.
The figures below illustrate the current vector potential, \(T\), the free-current density, \(\vec{J}_f\), and magnetic field, \(B\), computed with the following settings: \(p = 2\) and \(r = 4\). Visual inspection of the magnetic vector potential, \(\vec{A}\), is not very informative as its conservative portion is unknown.
The images above suggest that the computed fields do not exhibit any irregular behavior (the first image in the table below illustrates how irregular behavior can look like). The fields on these images closely resemble the corresponding closed-form analytical expressions given in the introduction. From the first glans the convergence tables above may appear somewhat strange due to extra rapid convergence at the higher degrees of the finite elements. One, however, can argue that the extra rapid convergence can be explained by the simple polynomial form of the fields being approximated. One thing is certain - the convergence rates presented in these tables are at the best theoretically expected values or better.
Let us consider the two-dimensional curl-curl partial differential equation again,
\begin{equation} \vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg) + \eta^2 \vec{A} = \vec{\nabla}\overset{V}{\times} T. \end{equation}
We can compute the divergence of this expression and rearrange the terms as the following:
\begin{equation} \vec{\nabla}\cdot \vec{A} = \frac{1}{\eta^2} \vec{\nabla}\cdot \bigg[\vec{\nabla}\overset{V}{\times} T - \vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg)\bigg]. \end{equation}
The divergence of the vector curl equals zero, see the introduction. Therefore, the right-hand side of the last equation equals zero,
\begin{equation} \vec{\nabla}\cdot \vec{A} = 0. \end{equation}
That is to say, the \(\eta^2\) gauging term can be considered to be the Coulomb gauge at least in theory. In practice the situation is a bit more complicated. What the \(\eta^2\) gauging term does depends on the value of \(\eta^2\). The table below attempts to express this very point.
This table presents four simulations for four different values of \(\eta^2\). In all four simulations the degree of the finite elements and the mesh refinement parameter were \(p=2\) and \(r=1\), respectively. The setting Settings::project_exact_solution was set to true.
At the setting \(\eta^2=\dfrac{10^{-10}}{\mu_0}\) the computed magnetic vector potential (the orange curve on the \(|\vec{A}(x,0)|\) plot) looks exactly the same as the projected exact solution gauged by the Coulomb gauge (the blue curve). At this value of \(\eta^2\) the term \(\eta^2\vec{A}\) acts like the Coulomb gauge. At the default setting, \(\eta^2=0\), the solution, \(\vec{A}\), is contaminated by an unknown conservative vector field. The computed magnetic field, \(B\), at this setting looks exactly as the corresponding exact solution because the conservative portion of \(\vec{A}\) is filtered out by the process of computing the magnetic field, \(B = \vec{\nabla}\overset{S}{\times}\vec{A}\). At the setting \(\eta^2=\dfrac{10^{-6}}{\mu_0}\) the solenoidal part of the computed magnetic vector potential deviates from the exact solution. We can deduce this from the fact that the computed magnetic field deviates from the exact expression for the magnetic field on the \(B(x,0)\) plot. The error in the solenoidal part of the computed \(\vec{A}\) is due to the fact that introduction of the gauging term, \(\eta^2\vec{A}\), modifies the initial curl-curl equation. So, strictly speaking, we are solving a different partial differential equation. If \(\eta^2\) is small, the difference between the initial curl-curl equation and the curl-curl equation modified by adding the gauging term is negligible. Consequently, the error in the solenoidal part of the computed \(\vec{A}\) is negligible as well. Evidently, in this particular case "small" means \(\eta^2 \ll \dfrac{10^{-6}}{\mu_0}\).
Normally, we are interested in measurable fields such as magnetic field and consider the magnetic vector potential as a useful tool for computing measurable fields. We can tolerate the presence of an unknown conservative component in the magnetic vector potential, i.e., the situation illustrated by the first two rows in the table above. In such disposition we need to keep \(\eta^2\) as small as possible, i.e., as far away as possible from the situation shown in last row of the table. In this tutorial program we set \(\eta^2\) to zero and increase it just a bit in if the conjugate gradient algorithm cannot converge.
Suppose for a moment that we would like to have the solution to the curl-curl equation in terms of a purely solenoidal magnetic vector potential, \(\vec{A}\), that is, a solution gauged by the Coulomb gauge. To get such a solution we need to tweak the \(\eta^2\) parameter. By contemplating the table above one can hypothesize that there exists an optimal value of the gauging parameter, \(\eta^2_\text{opt}\), at which the \(L^2\) error norm computed for \(\vec{A}\) is minimal. The optimal value should be somewhere in between \(\eta^2=\dfrac{10^{-12}}{\mu_0}\) and \(\eta^2=\dfrac{10^{-6}}{\mu_0}\), based on the experiments above. Try to verify this hypothesis by finding the exact value of \(\eta^2_\text{opt}\).
In this instance we have a close-form analytical expression of the exact solution, i.e., the expression of \(\vec{A}\) given in the introduction. In a real-life simulation there is no expression of the exact solution. Try to think of a method of blind (meaning without the exact solution) estimation of \(\eta^2_\text{opt}\). Try to implement and test your ideas.
Adding the gauging term, \(\eta^2\vec{A}\), to the curl-curl equation converts a positive semidefinite system matrix into a positive definite matrix. Strictly speaking, adding the gauging term modifies the initial curl-curl equation. For this reason, it is important to keep \(\eta^2\) small so the solution is not afflicted by the error induced by adding the gauging term. The definition of "small" here is a bit fuzzy. In absence of the exact solution setting \(\eta^2\) at the acceptable level of the error in the solenoidal component of \(\vec{A}\) is difficult. It is better to discard the gauging term and use a more sophisticated linear solver. The hypre AMS [137] can solve the systems of linear equations yielded by the curl-curl equation without the \(\eta^2\) gauging term, ( \(\beta=0\) in [137]). Try to implement the hypre AMS.
1.9.8