Introduction

In this numerical experiment we will calculate the magnetic field, $\vec{B}$, induced by a thick spherical coil. The figure below illustrates a cross section of the coil.

The calculation is done in four stages:

Calculate the current vector potential, $\vec{T}$, for a given volume free-current density $\vec{J}_f$ (StaticVectorSolver::Solver1).
Check the calculated $\vec{T}$ by converting it back into $\vec{J}_f$ (StaticVectorSolver::ProjectTtoJ).
Use the calculated current vector potential to solve for $\vec{A}$ (StaticVectorSolver::Solver2).
Convert computed $\vec{A}$ into the magnetic field $\vec{B}$ (StaticVectorSolver::ProjectAtoB).

The Bossavit's diagrams in the figure below illustrate these four stages.

This numerical experiment is carried out in the three-dimensional space. As implied by the Bossavit's diagrams above, see also this discussion, the vector potentials, $\vec{T}$ and $\vec{A}$, belong to the $H(\text{curl})$ function space, and, thus, must be modeled by the FE_Nedelec finite elements. Similarly, $\vec{J}_f$ and $\vec{B}$ belong to the $H(\text{div})$ function space and, thus, must be modeled by the FE_RaviartThomas finite elements.

As per usual, we will execute the experiment on a family of progressively refined meshes and compare the numerical results to closed-form analytical expression of the exact solution. We will vary the degree of the finite elements and observe that the higher degree of the finite elements imply faster convergence of the $L^2$ error norms. As soon as we do not gauge the vector potentials ( $\vec{T}$ and $\vec{A}$) explicitly, their conservative parts will remain unknown. For this reason, it does not make much sense to visualize them or evaluate them by computing error norms. We will evaluate the computed vector potentials indirectly by computing $L^2$ error norms of the fields derived from them ( $\vec{J}_f$ and $\vec{B}$).

Implementation

The problem and the problem domain

The following textbook problem is considered (Section 6.2.5. Thick spherical coil).

Volume free-current density, $\vec{J}_f$, is confined to an interior of a spherical shell. The shell is shown in blue color in the figure above. The inner and outer radii of the shell equal $ a $ and $ b $, respectively. The free-current density can be expressed in a spherical coordinate system as
\[ \vec{J}_f = \left\{ \begin{aligned} & 0 &&\text{if }&& r < a\\ & K_0 r \sin(\theta) \hat{\phi} &&\text{if }&& a \le r \le b\\ & 0 &&\text{if }&& r > b.\\ \end{aligned} \right. \]

We would like to calculate the magnetic field, $\vec{B}$, induced by the coil observing the compatibility condition for the curl-curl equation.

The figure below illustrates the geometry of the coil and the corresponding three-dimensional problem domain.

The Dirichlet boundary condition,

\[ \hat{n} \times \vec{T} = 0, \]

is applied when computing the current vector potential. The boundary ID in this case equals $1$. This boundary ID is set in the geo files. The Robin boundary condition (first-order ABC),

\[ \frac{1}{\mu} \hat{n} \times \bigg(\vec{\nabla} \times \vec{A}\bigg) + \frac{1}{\mu r} \hat{n} \times \bigg(\hat{n} \times \vec{A}\bigg) = 0, \]

is applied when computing the magnetic vector potential. The boundary ID in this case equals $2$. This boundary ID is set by the code of the program. That is, during the simulation the boundary ID has two values: $1$ and $2$.

The exact solution

The problem has the following closed-form analytical solution.

\[ \vec{B} = \left\{ \begin{aligned} & \frac{1}{2} (b^2-a^2) \vec{F}_1(\theta) &&\text{if }&& r \le a\\ & \frac{1}{2} (b^2-r^2) \vec{F}_1(\theta) + \frac{1}{5} (r^5-a^5) \vec{F}_2(r, \theta) && \text{if } && a \le r \le b\\ & \frac{1}{5} (b^5-a^5) \vec{F}_2(r, \theta) &&\text{if }&& r \ge b,\\ \end{aligned} \right. \]

where

\[ \begin{aligned} &\vec{F}_1(\theta) = \dfrac{2}{3}\mu_0 K_0\big(\cos(\theta)\hat{r}-\sin(\theta)\hat{\theta}\big), \\ &\vec{F}_2(r, \theta) = \dfrac{2}{3}\mu_0 K_0 \frac{1}{r^3}\big(\cos(\theta)\hat{r}+\frac{1}{2}\sin(\theta)\hat{\theta}\big). \end{aligned} \]

The mesh

The figure below illustrates the mesh.

The main elements of the mesh are listed below.

The cube with a side of $2d_1$ in the middle of the mesh. The purpose of this cube is instrumental. This is the only way to construct a mesh that has spherical interfaces and contains no tetrahedral cells as the cells around the origin must be hexahedral.
The spherical shell. This region contains the free current, $ \vec{J}_f = K_0 r \sin(\theta) \hat{\phi} $. There is no free current outside this region.
The sphere nr.1. The $L^2$ error norms are computed inside the region delineated by this surface.
The sphere nr.2. The boundary of the problem domain. Represents infinity. The boundary conditions (Dirichlet, Robin) is applied to this surface.

The boundary value problem

At the 0th stage of the experiment the current vector potential, $\vec{T}$, is calculated by solving the following boundary value problem

\begin{equation} \begin{array}{lrcll} \text{ } & \vec{\nabla}\times\bigg(\vec{\nabla}\times\vec{T}\bigg) +\eta^2\vec{T} = \vec{\nabla}\times\vec{J}_f & \text{in} & \Omega & \text{(i)},\\ \text{(e)}& \hat{n} \times \vec{T} = 0 & \text{on} & \Gamma_{R1} & \text{(ii)}. \\ \end{array} \end{equation}
At the 1st stage of the experiment the current vector potential calculated at the 0th stage, $\vec{T}$, is converted back into free-current density as
\[ \vec{J}_{f} = \vec{\nabla} \times \vec{T} \]
for the purpose of evaluation.
At the 2nd stage of the experiment the following boundary value problem is solved:
\begin{equation} \begin{array}{lrcll} \text{ } & \vec{\nabla}\times\bigg(\dfrac{1}{\mu_0} \vec{\nabla}\times\vec{A}\bigg) +\eta^2\vec{A} = \vec{\nabla} \times \vec{T} & \text{in} & \Omega & \text{(i)},\\ \text{(n)}& \dfrac{1}{\mu_0} \hat{n} \times \bigg(\vec{\nabla} \times \vec{A}\bigg) + \gamma \hat{n} \times \bigg(\hat{n} \times \vec{A}\bigg) = 0 & \text{on} & \Gamma_{R1} & \text{(ii)},\\ \end{array} \end{equation}
where $ \vec{T} $ is the current vector potential computed at the 0th stage and
\[ \gamma = \frac{1}{\mu_0 r}. \]
At the 3rd stage of the experiment the magnetic vector potential calculated at the 2nd stage, $\vec{A}$, is converted into magnetic field as:
\[ \vec{B} = \vec{\nabla} \times \vec{A}. \]

The program

This numerical experiment is implemented in accordance with the base structure described in here. All computations are done in four stages. The stages are labeled from 0 to 3. The following algorithms where used:

At the 0th stage the current vector potential, $\vec{T}$, is computed with a help of the StaticVectorSolver::Solver1 class template.
At the 1st stage the current vector potential computed at the 0th stage is converted back into volume free-current density, $\vec{J}_{f}$, with a help of the StaticVectorSolver::ProjectTtoJ class template.
At the 2nd stage the magnetic vector potential, $\vec{A}$, is computed with a help of the StaticVectorSolver::Solver2 class template. The current vector potential, $\vec{T}$, computed at the 0th stage is fed to the algorithm.
At the 3rd stage the magnetic vector potential computed at the 2nd stage is converted to the magnetic field, $\vec{B}$, with a help of the StaticVectorSolver::ProjectAtoB class template.

The Dirichlet and Robin (ABC) boundary conditions are applied at the 0th and the 2nd stages, respectively. The boundary ID must satisfy the boundary ID convension. For this reason, two boundary ID's are used. At the 0th stage the boundary ID equals $1$. This boundary ID is set in the geo file as the following:

Physical Surface(1) = {9, 189, 433, 309, 553, 673};

It is loaded together with the mesh. The boundary ID equals $2$ at the second stage. The boundary ID is changed from $1$ to $2$ in the constructor of the SolverSSOLII_A class:

...
 
SolverSSOLII_A(unsigned int p,
               unsigned int mapping_degree,
               unsigned int r,
               const Triangulation<3>& triangulation_T,
               const DoFHandler<3>& dof_handler_T,
               const Vector<double>& solution_T,
               std::string fname)
  : Solver2<3, 2>(p,
                  mapping_degree,
                  triangulation_T,
                  dof_handler_T,
                  solution_T,
                  2,
                  0.0,
                  fname,
                  nullptr,
                  SettingsSSOLII::print_time_tables,
                  false,
                  true)
  , r(r)
  , fname(fname)
{
  for (auto cell : Solver2<3, 2>::triangulation_T.active_cell_iterators()) {
    for (unsigned int f = 0; f < GeometryInfo<3>::faces_per_cell; f++) {
      if (cell->face(f)->at_boundary() && (cell->face(f)->boundary_id() == 1))
        cell->face(f)->set_boundary_id(2);
    }
  }
 
  Solver2<3, 2>::run();
}
 
...

Note that the free-current density, $\vec{J}_f$, equals zero on the boundary by definition of the problem. The tangential component of the current vector potential, $\vec{T}$, equals zero on the boundary as we apply the homogeneous Dirichlet boundary condition when computing it, see the first boundary value problem above. For these two reasons we can neglect integrals $I_{b3-2}$ in the Recipe for static vector solver in 3D (current vect. potential) and in Recipe for static vector solver in 3D. To do so, we pass type_of_pde_rhs = 2 to the constructors of both, StaticVectorSolver::Solver1 and StaticVectorSolver::Solver2.

The build process generates one executable file: ssol-ii-sphere. To rebuild it change into ssol-ii/build/Release directory and execute the following:

user@computer .../ssol-ii/build/Release$ ./clean

user@computer .../ssol-ii/build/Release$ ./build

Then the executable file must be executed again. This can be done as the following:

user@computer .../ssol-ii/build/Release$ cd ../../bin/Release

user@computer .../ssol-ii/bin/Release$ ./ssol-ii-sphere

This will generate vtu files in the ssol-ii/bin/Release/Data directory. They can be viewed with a help of a fresh version of ParaView. The VisIt visualization software will not do as this numerical experiment uses second-degree mapping.

Note that executable files require mesh files to be present in the ssol-ii/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into ssol-ii/gmsh directory and executing the following:

user@computer .../ssol-ii/gmsh$ ./clean

user@computer .../ssol-ii/gmsh$ ./build

This will generate a set of globally refined meshes in ssol-ii/gmsh/data directory.

The SettingsSSOLII class allows switching on three useful features: printing time tables on the computer screen, logging convergence data of the conjugate gradient solver, and saving the exact solution into the vtu files next to the numerical solution. The last can be useful when debugging.

Simulation results

The simulation was done with the following settings: $d_1 = 0.1[m]$, $a = 0.9[m]$, $b = 1.2[m]$, $d_2 = 2.0[m]$, $d_3 = 3.0[m]$, and $K_0=1[A\cdot m^{-2}]$.

The first figure on this page illustrates a vector plot of the computed free-current density, $\vec{J}_f$, and the magnetic field, $\vec{B}$.

The free-current density

The figure below illustrates two slices of the magnitude of free-current density, $|\vec{J}_{f}|$, calculated with finite elements of degree $p = 2$ on the mesh sphere_r9.msh.

The table below illustrates the convergence of the $L^2$ norm of the error made in computing the free-current density, $\vec{J}_{f}$.

Free-current density, $\vec{J}_{f}$

p	r	cells	dofs	$\\|e\\|_{L^2}$	$\alpha_{L^2}$
0	6	3125	9450	1.66e-01	-
0	7	5400	16308	1.38e-01	0.99
0	8	8575	25872	1.19e-01	0.99
0	9	12800	38592	1.04e-01	0.99
1	6	3125	75300	8.12e-04	-
1	7	5400	130032	4.97e-04	2.69
1	8	8575	206388	3.32e-04	2.61
1	9	12800	307968	2.37e-04	2.54
2	6	3125	253800	6.78e-04	-
2	7	5400	438372	3.94e-04	2.97
2	8	8575	695898	2.49e-04	2.98
2	9	12800	1038528	1.67e-04	2.99

The following notations were used in the table header:

p - the degree of the FE_Nedelec finite elements that were used to model the current vector potential, $\vec{T}$, and the degree of the FE_RaviartThomas finite elements that were used to model the free-current density, $\vec{J}_f$.
r - the number of nodes on transfinite lines, see discussion in here.
cells - the total amount of active cells.
dofs - the amount of degrees of freedom in the active cells.
$\|e\|_{L^2}$ - the $L^2$ error norm.
$\alpha_{L^2}$ - the order of convergence of the $L^2$ error norm.

The magnetic field

The figure below illustrates two slices of the magnitude of magnetic field, $|\vec{B}|$, calculated with finite elements of degree $p = 2$ on the mesh sphere_r9.msh.

The table below illustrates the convergence of the $L^2$ norm of the error made in computing the magnetic field, $\vec{B}$.

p	r	cells	dofs	$\\|e\\|_{L^2}$	$\alpha_{L^2}$
0	6	3125	9450	7.12e-08	-
0	7	5400	16308	5.91e-08	1.02
0	8	8575	25872	5.06e-08	1.01
0	9	12800	38592	4.40e-08	1.04
1	6	3125	75300	4.75e-09	-
1	7	5400	130032	3.68e-09	1.40
1	8	8575	206388	2.85e-09	1.65
1	9	12800	307968	2.41e-09	1.27
2	6	3125	253800	1.23e-10	-
2	7	5400	438372	6.53e-11	3.46
2	8	8575	695898	3.57e-11	3.92
2	9	12800	1038528	2.27e-11	3.39