Introduction

The purpose of this numerical experiment is to test the StaticVectorSolver::Solver2 class template in the two-dimensional space by the method of manufactured solutions. The class templates StaticVectorSolver::Solver1 and StaticVectorSolver::Solver2 are very similar. Both class templates solve the static vector boundary value problem. In two spatial dimensions (dim=2) these two solvers differ from each other only at two points:

When solving a problem formulated in terms of the magnetic vector potential, $\vec{A}$, the algorithm implemented by StaticVectorSolver::Solver1 class template assumes that the right-hand side of the partial differential equation is either volume free-current density, $\vec{J}_f$, or the curl of the current vector potential, $\vec{\nabla}\overset{V}{\times} T$. In contrast, the StaticVectorSolver::Solver2 class template allows only $\vec{\nabla}\overset{V}{\times}T$ on the right-hand side of the partial differential equation.
In the StaticVectorSolver::Solver1 class template the field on the right-hand side of the partial differential equation, i.e., $\vec{J}_f$ or $T$, must be given as a closed-form analytical expression. In contrast, the StaticVectorSolver::Solver2 class template assumes that the vector field on the right-hand side of the partial differential equation, $T$, is a field function, i.e., is a result of another deal.II simulation.

All other aspects of the class templates StaticVectorSolver::Solver1 and StaticVectorSolver::Solver2 are the same. Due to this similarity most of the aspects of this numerical experiment are exactly the same as that of the two-dimensional version of the Method of manufactured solutions, vector potential (mms-v/) numerical experiment. The difference is in how the vector field on the right-hand side of the partial differential equation, $T$, is fed to the solver. In the Method of manufactured solutions, vector potential (mms-v/) numerical experiment $T$ is fed to the solver in a form of analytical expression. In this numerical experiment $T$ is computed numerically in a preprocessing step and then fed to the solver that computes $\vec{A}$.

Another difference between the Method of manufactured solutions, vector potential (mms-v/) numerical experiment and this numerical experiment is in the amount of spatial dimensions in the solved problems. In the Method of manufactured solutions, vector potential (mms-v/) numerical experiment both, two- and three-dimensional problems are considered. In this numerical experiment - only two-dimensional. Test of the StaticVectorSolver::Solver2 class template in the three-dimensional space is done in the Method of manufactured solutions, vector potential (mms-vt-i) numerical experiment.

In this numerical experiment we first manufacture a solution, $\vec{A}_M$. After that, we differentiate it to get a closed-form analytical expression for the surface free-current density, $\vec{J}_M$, and the vector fields on the right-hand sides of the static vector boundary value problem. We choose the coefficients on the left-hand sides somewhat arbitrary. Next, to solve for $\vec{A}$ we do the following:

Calculate numerically the current vector potential, $T$, for a given $\vec{J}_M$ (the preprocessing step).
Invoke the algorithm of StaticVectorSolver::Solver2 to solve for $\vec{A}$. Use $T$ computed in the preceding step as an input.
Convert computed $\vec{A}$ into the magnetic field $B$ (the postprocessing step).

The last step, i.e., converting $\vec{A}$ into the magnetic field $B$, is necessary as we do not gauge the magnetic vector potential explicitly. The last implies that the magnetic vector potential is not unique and we do not know in advance which of the many possible vector potentials will appear at the output of the solver. Consequently, we cannot compare the calculated magnetic vector potential with the manufactured magnetic vector potential. Instead we will compare the calculated and manufactured magnetic fields. The numerically calculated magnetic field is unique even if the magnetic vector potential is ungauged.

Having computed $B$ we compare it to the manufactured magnetic field,

\[ B_M = \vec{\nabla} \overset{S} \times \vec{A}_M, \]

by computing an error as

\[ e = | B_M-B | \]

and evaluating its $L^2$ error norm as

\[ ||e||_{L^2} = \sqrt{\iint_{\Omega} e^2 dS}. \]

The three simulation steps above, i.e., computing $T$, $\vec{A}$, and $B$, are mandatory. We add to them one optional simulation for checking if the numerically calculated current vector potential, $T$, has been calculated correctly. To this end we convert the numerically calculated current vector potential back to the volume free-current density,

\[ \vec{J}_{f} = \vec{\nabla} \overset{V} \times T, \]

compute an error,

\[ e = |\vec{J}_{M} - \vec{J}_{f}|, \]

where

\[ \vec{J}_M = \vec{\nabla}\overset{V}\times\bigg(\frac{1}{\mu}\vec{\nabla}\overset{S} \times\vec{A}_M\bigg) \]

is the manufactured volume free-current density, and then evaluate the $L^2$ error norm,

\[ ||e||_{L^2} = \sqrt{\iint_{\Omega} e^2 dS}. \]

Note that in two-dimensional problems the current vector potential is implicitly gauged by the translation symmetry. For this reason we can directly compare the numerically computed current vector potential, $T$, to its manufactured counterpart, $T_M$, by computing the $L^2$ and $H^1$ error norms.

The four computation steps (stages) discussed above are illustrated by means of the Bossavit's diagrams in the figure below.

Implementation

The manufactured solution and the problem domain

We choose (somewhat arbitrary) the following manufactured solution:

\[ \vec{A}_M = \frac{1}{k}\big(-\sin{(ky)}\hat{i} + \sin{(kx)} \hat{j}\big). \]

Note that

\[ \vec{\nabla} \cdot \vec{A}_M = 0. \]

Next, we calculate the manufactured magnetic field as

\[ B_M = \vec{\nabla} \overset{S}{\times} \vec{A}_M = \begin{vmatrix} 0 & 0 & 1 \\ \dfrac{\partial}{\partial x}&\dfrac{\partial}{\partial y}& 0 \\ -\dfrac{1}{k}\sin{(ky)} & \dfrac{1}{k}\sin{(kx)} & 0 \end{vmatrix} = \cos{(kx)} + \cos{(ky)}. \]

We choose the coefficients on the left-hand side of the boundary value problem somewhat arbitrary as well,

\[ \mu = \mu_0 \big( x^2 + y^2 + 1 \big), \]

\[ \gamma = \frac{1}{\mu} \bigg( \sqrt{x^2 + y^2} + 2 \bigg). \]

Next, we calculate the manufactured free-current density,

\[ \begin{aligned} & \vec{J}_M = \vec{\nabla}\overset{V}\times\bigg(\frac{1}{\mu}\vec{\nabla}\overset{S}\times\vec{A}_M\bigg) = \vec{\nabla}\overset{V}\times\bigg(\frac{1}{\mu}B_M\bigg) = \frac{\partial}{\partial y} \bigg(\frac{1}{\mu} B_M \bigg) \hat{i} - \frac{\partial}{\partial x} \bigg(\frac{1}{\mu} B_M \bigg) \hat{j} = \\ & = \bigg( B_M \frac{\partial}{\partial y} \frac{1}{\mu} + \frac{1}{\mu} \frac{\partial}{\partial y} B_M \bigg) \hat{i} - \bigg( B_M \frac{\partial}{\partial x}\frac{1}{\mu} + \frac{1}{\mu}\frac{\partial}{\partial x} B_M \bigg) \hat{j}. \end{aligned} \]

Next, we note that

\[ \frac{\partial}{\partial x} \frac{1}{\mu} = \frac{1}{\mu_0} \frac{\partial}{\partial x} (x^2+y^2+1)^{-1} = -\frac{1}{\mu_0} (x^2+y^2+1)^{-2} (2x) = -\frac{2x}{\mu_0 (x^2+y^2+1)^2} = - 2 x \frac{\mu_0}{\mu^2} \]

and

\[ \frac{\partial}{\partial y} \frac{1}{\mu} = - 2 y \frac{\mu_0}{\mu^2}. \]

Therefore, we can rewrite the manufactured volume free-current density as

\[ \begin{aligned} & \vec{J}_M = \bigg(-B_M 2 y \frac{\mu_0}{\mu^2} + \frac{1}{\mu} \frac{\partial}{\partial y} B_M \bigg) \hat{i} - \bigg(-B_M 2 x \frac{\mu_0}{\mu^2} + \frac{1}{\mu} \frac{\partial}{\partial x} B_M \bigg) \hat{j} = \\ &=\frac{1}{\mu}\Bigg[ \bigg(- 2 y \frac{\mu_0}{\mu}\big(\cos{(kx)} + \cos{(ky)}\big) - k \sin{(ky)} \bigg) \hat{i} - \bigg(- 2 x \frac{\mu_0}{\mu}\big(\cos{(kx)} + \cos{(ky)}\big) - k \sin{(kx)} \bigg) \hat{j} \Bigg]. \end{aligned} \]

The last equation is used for calculating the $L^2$ error norms as discussed in the introduction.

Finally, we calculate the right-hand sides of the boundary value problem as

\[ \vec{G} = \vec{A}_M \]

and

\[ \vec{Q} = \frac{1}{\mu}\hat{n}\overset{V}{\times}\bigg(\vec{\nabla}\overset{S}{\times}\vec{A}_M\bigg) +\gamma \hat{n} \overset{V}{\times} \bigg(\hat{n}\overset{S}{\times}\vec{A}_M\bigg) = \hat{n} \overset{V}{\times} T_M + \gamma \hat{n} \overset{V}{\times} \bigg(\hat{n}\overset{S}{\times}\vec{A}_M\bigg), \]

where

\[ T_M = \frac{1}{\mu} \vec{\nabla} \overset{S}{\times} \vec{A}_M = \frac{1}{\mu} B_M. \]

Two problem domains are used. The first problem domain is shaped as a square The second problem domain is shaped as a circle. In both problem domains the Dirichlet boundary condition is applied to a half of the boundary. The Robin boundary condition is applied to another half. The boundary IDs are assigned in accordance with the boundary ID convention. The two problem domains are shown in the figure below.

The boundary value problem

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

\begin{equation} \begin{array}{lrcll} \text{ } & - \vec{\nabla}\cdot\big(\vec{\nabla} T \big) = \vec{\nabla}\overset{S}{\times}\vec{J}_M & \text{in} & \Omega & \text{(i)},\\ \text{(e)}& T = T_M & \text{on} & \Gamma_{D1} & \text{(ii)}, \\ \text{(n)}& \hat{n} \cdot \vec{\nabla} T = -n_x J_{My} + n_y J_{Mx} & \text{on} & \Gamma_{R1} & \text{(iii)}, \end{array} \end{equation}

where $ T_M $ is the manufactured current vector potential and $J_M$ is the manufactured free-current density. Note that
\[ \vec{J}_f = \vec{\nabla}\overset{V}{\times}T = \frac{\partial T}{\partial y} \hat{i} - \frac{\partial T}{\partial x} \hat{j}. \]

Therefore,
\[ \hat{n}\cdot\vec{\nabla} T = n_x\frac{\partial T}{\partial x} + n_y\frac{\partial T}{\partial y}= -n_x J_{fy} + n_y J_{fx}. \]

Thus the equation (iii) in the boundary value problem above.
At the 1st stage the current vector potential computed at the 0-th stage, $ T $, is converted into the volume free-current density as
\[ \vec{J}_{f} = \vec{\nabla}\overset{V}{\times} T. \]
Then computed $\vec{J}_{f}$ is compared to the manufactured free-current density $\vec{J}_M$.
At the 2nd stage the magnetic vector potential $\vec{A}$ is computed by solving the following boundary value problem

\begin{equation} \begin{array}{lrcll} \text{ } & \vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg) = \vec{\nabla} \overset{V}{\times} T & \text{in} & \Omega & \text{(i)},\\ \text{(e)}&\hat{n}\overset{S}{\times}\vec{A} = \hat{n}\overset{S}{\times}\vec{G} & \text{on} & \Gamma_{D1} & \text{(ii)}, \\ \text{(n)} &\dfrac{1}{\mu} \hat{n} \overset{V} \times \bigg(\vec{\nabla}\overset{S}\times \vec{A} \bigg) + \gamma \hat{n}\overset{V}{\times} \bigg( \hat{n}\overset{S}{\times}\vec{A} \bigg) = \vec{Q} & \text{on} & \Gamma_{R1} & \text{(iii)}, \\ \end{array} \end{equation}

where $ T $ is the current vector potential computed at the 0-th stage.
At the 3rd stage the magnetic vector potential computed at the 2nd stage, $\vec{A}$, is converted into magnetic field as
\[ B = \vec{\nabla}\overset{S}{\times} \vec{A}. \]
Then $ B $ is compared to the manufactured magnetic field, $ B_M $, by computing the $L^2$ error norm, see the introduction.

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, see the preceding section. The following class templates where used:

At the 0th stage the current vector potential, $T$, is computed with a help of the StaticScalarSolver::Solver 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 StaticScalarSolver::ProjectTzToJxy 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, $T$, computed in the 0th stage is fed to the algorithm as the scalar field under the vector curl on the right-hand side of the partial differential equation.
At the 3rd stage the magnetic vector potential computed in the 2nd stage is converted to the magnetic field, $B$, with a help of the StaticVectorSolver::ProjectAxyToBz class template.

Note that unlike most of the numerical experiments, this numerical experiment has two input files: mms-vt-ii/src/static_vector_input.cpp and mms-vt-ii/src/static_scalar_input.cpp.

The build process generates two executable files: mms-vt-ii-square, mms-vt-ii-circle. To rebuild them change into mms-vt-ii/build/Release directory and execute the following:

user@computer .../mms-vt-ii/build/Release$ ./clean

user@computer .../mms-vt-ii/build/Release$ ./build

Then all executable files must be executed again. This can be done by changing into the mms-vt-ii/bin/Release directory and executing run-all script there,

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

user@computer .../mms-vt-ii/bin/Release$ ./run-all

This will generate vtu files in the mms-vt-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. The data files in tex and txt format contain the convergence tables.

Note that executable files require a set of meshes to be present in the mms-vt-ii/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into mms-v/gmsh directory and executing the following:

user@computer .../mms-vt-ii/gmsh$ ./clean

user@computer .../mms-vt-ii/gmsh$ ./build

This will generate a set of refined meshes in mms-vt-ii/gmsh/data.

The SettingsMMSVTII 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 (manufactured) solutions into the vtu files next to the numerical solution. The last can be useful when debugging.

Simulation results

All the simulations described below have been done with a setting $k=\pi$. Recall, we do not gauge the magnetic vector potential, $\vec{A}$, explicitly. For this reason, presenting simulation results in terms of $\vec{A}$ does not make any sense: the solution being expressed in terms of $\vec{A}$ is not unique. Instead, we present the results in terms of the curl of the magnetic vector potential, i.e., in terms of $\vec{B}$.

Square domain

The three figures below present the results of the simulations on the square domain.

The corresponding convergence tables are presented below.

The following notations were used in the table headers:

p - the degree of the FE_RaviartThomas or FE_DGQ finite elements. The degree of the FE_Q finite elements is $p+1$
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.
$\|e\|_{H^1}$ - the $H^1$ error norm.
$\alpha_{L^2}$ - the order of convergence of the $L^2$ error norm.
$\alpha_{H^1}$ - the order of convergence of the $H^1$ error norm.

Circular domain

The three figures below present the results of the simulations on the circular domain.