Introduction

Some problems in magnetostatics can be formulated such that all free currents, $\vec{J}_f$, are left outside the problem domain. In such cases the auxiliary vector field $\vec{H}$ can be derived as a negative gradient of the total magnetic scalar potential,

\[ \vec{H} = - \vec{\nabla} \Psi \]

(see Section 3.1. Total magnetic scalar potential). In this case the behavior of the auxiliary vector field $\vec{H}$ is similar to the behavior of the electrostatic field $\vec{E}$. This similarity is not surprising as the electrostatic field can be derived from the electrostatic scalar potential in the exact same way and is governed by the exact same partial differential equation.

This numerical experiment is made as an illustration of application of the total magnetic scalar potential, $\Psi$. Two problems are considered: two-dimensional and three-dimensional. We will solve these two problems numerically on a set of progressively refined meshes and compare the numerical solutions to closed-form solutions derived by analytical methods. As per usual, we hope that the convergence rates will decrease as discussed in here

Implementation

Two-dimensional problem

The following textbook problem is used in the two-dimensional version of the experiment (Section 6.2.3. Cylindrical shield in a uniform magnetic field).

A magnetostatic shield is made of soft paramagnetic material of a permeability $\mu_1$. The magnetic material is lossless, linear, homogeneous, and isotropic. The shield is shaped as an infinitely long cylindrical tube coaxial with the $z$ axis. The inner and outer radii of the shield are $a$ and $b$, respectively. The space inside and outside the shield is free, so the permeability of the entire space can be described as
\begin{equation} \mu = \left\{ \begin{aligned} & \mu_0 &\text{if } \text{ } & r < a_1\\ & \mu_1 &\text{if } \text{ } & a_1 \le r \le b_1\\ & \mu_0 &\text{if } \text{ } & r > b_1, \end{aligned} \right. \end{equation}

where $r = \sqrt{x^2+y^2}$. The shield is exposed to a uniform magnetic field. The magnetic field in absence of the shield is described by the following equation
\[ \vec{B} = B_0 \hat{i}. \]

The figure below illustrates this configuration and the corresponding problem domain.

We would like to calculate the total magnetic scalar potential, $\Psi$, the magnetic field, $\vec{B}$, and the auxiliary vector field $\vec{H}$ in a two-dimensional planar problem domain.

The closed-form solution to this problem obtained by analytical methods reads

\begin{equation} \Psi = \left\{ \begin{aligned} &\delta_1 x &\text{if } \text{ } & r \le a\\ &\bigg(\beta_1+\frac{\gamma_1}{r^2}\bigg) x &\text{if } \text{ } & a \le r \le b\\ &\bigg(-H_0 + \frac{\alpha_1}{r^2}\bigg) x &\text{if } \text{ } & r \ge b, \end{aligned} \right. \end{equation}

\begin{equation} \vec{\nabla} \Psi = \left\{ \begin{aligned} &\delta_1 \hat{i} &\text{if } \text{ } & r \le a\\ & \bigg(\beta_1 + \gamma_1\frac{1}{r^2} - 2\gamma_1 \frac{x^2}{r^4} \bigg)\hat{i} -2\gamma_1\frac{xy}{r^4}\hat{j} &\text{if } \text{ } & a \le r \le b\\ & \bigg(-H_0 + \alpha_1\frac{1}{r^2} - 2\alpha_1\frac{x^2}{r^4} \bigg)\hat{i} -2\alpha_1\frac{xy}{r^4}\hat{j} &\text{if } \text{ } & r \ge b, \end{aligned} \right. \end{equation}

\[ \vec{H} = - \vec{\nabla} \Psi, \]

\[ \vec{B} = - \mu \vec{\nabla} \Psi, \]

where

\[ r = \sqrt{x^2 + y^2}, \]

\[ \mu_r = \frac{\mu_1}{\mu_0}, \]

\begin{equation} \begin{aligned} &\Omega = \frac{(\mu_r - 1)}{(\mu_r + 1)} \frac{a^2}{b^2}, \\ &\gamma_1 = \frac{-2b^2 H_0 \Omega}{(\mu_r+1)-(\mu_r-1) \Omega}, \\ &\beta_1 = \frac{(\mu_r+1) \gamma_1}{(\mu-1)a^2}, \\ &\alpha_1 = -b^2 H_0+\mu_r\gamma_1-\mu_rb^2\beta_1, \\ &\delta_1 = \frac{\mu_r a^2\beta_1-\mu_r\gamma_1}{a^2}. \end{aligned} \end{equation}

The mesh of the problem domain is shown in the figure below.

The list below describes the main elements of the mesh.

The disk a radius $a$. It represents the inner space of the shield. The permeability of the material inside the disk equals $\mu_0$.
The circular ring with an inner and outer radii of $a$ and $b$, respectively. It represents the shield. The permeability in this region equals $\mu_1$.
The region between the circle of a radius $b$ and the outer square boundary. This region represents the space outside the shield. The permeability in this region equals $\mu_0$.
The square with a side of $ 2d_2 $. This square delineates the local mesh. The error norms are computed inside the local mesh. The error norms outside the local mesh are discarded.
The square with a side of $2d_1$ in the middle of the mesh. The purpose of this square is instrumental. This is the only way of constructing a mesh that has circular interfaces and contains no triangular cells: the cells around the origin must be quadrangles.
The square with a side of $2d_3$ represents the infinity.

The H-field induced by the magnetization of the magnetic material of the shield vanish at infinity. Therefore, the total H-field at infinity equals the applied H-field. Thus, the Dirichlet boundary condition

\[ \Psi = -H_0 x, \]

where $H_0 = \dfrac{1}{\mu_0} B_0$, can be applied on all segments of the boundary. This corresponds to the "DirichletOnly" mode of the program. On the other hand, the applied H-field runs parallel the $x$ axis. Therefore, its normal component on the two horizontal boundaries equals zero. Consequently, we can apply the Dirichlet boundary condition on the left and right boundaries and the homogeneous Neumann boundary condition,

\[ \mu \hat{n}\cdot\vec{\nabla} \Psi = 0, \]

on the top and bottom boundaries. The program in the "DirichletNeumann" mode does exactly this. The problem domain shown above is truncated. That is, the boundary of the problem domain is at a finite distance from the origin. The domain of the initial problem is unbounded. For this reason, application of the Dirichlet boundary condition or a combination of the Dirichlet and Neumann boundary conditions will lead to a truncation error. We would like to observe the convergence rate without the truncation error. In general, the truncation error can be driven to be insignificantly small by increasing the parameter $d_3$, see the figure above. This, however, will increase the simulation time. Instead, we will keep the size of the problem domain relatively small and use the exact boundary conditions. That is, we will apply the Dirichlet boundary condition to all segments of the boundary and use the known exact solution (instead of $ -H_0 x $) to compute the right-hand side $\eta$. This will allow us to switch off the truncation error and observe the convergence rates free of the truncation error. This is what program does if switched into the "Exact" mode. The user can switch the program between the three modes by setting Setting::type_of_bc.

Three-dimensional problem

The following textbook problem is used in the three-dimensional version of the experiment. (Section 6.2.2. Spherical shield in a uniform magnetic field)

A magnetostatic shield is made of soft paramagnetic material of a permeability $\mu_1$. The magnetic material is lossless, linear, homogeneous, and isotropic. The shield is shaped as a spherical shell concentric with the origin. The inner and outer radii of the shield are $a$ and $b$, respectively. The space inside and outside the shield is free, so the permeability of the entire space can be described as
\begin{equation} \mu = \left\{ \begin{aligned} & \mu_0 &\text{if } \text{ } & r < a_1\\ & \mu_1 &\text{if } \text{ } & a_1 \le r \le b_1\\ & \mu_0 &\text{if } \text{ } & r > b_1, \end{aligned} \right. \end{equation}

where $r = \sqrt{x^2+y^2+z^2}$. The shield is exposed to a uniform magnetic field. The magnetic field in absence of the shield is described by the following equation
\[ \vec{B} = B_0 \hat{k}. \]

The figure below illustrates this configuration and the corresponding problem domain.

We would like to calculate the total magnetic scalar potential, $\Psi$, the magnetic field, $\vec{B}$, and the auxiliary vector field $\vec{H}$ in a three-dimensional problem domain.

The closed-form solution to this problem obtained by analytical methods reads

\begin{equation} \Psi = \left\{ \begin{aligned} &\delta_1 z &\text{if } \text{ } & r \le a\\ &\bigg(\beta_1 + \frac{\gamma_1}{r^3}\bigg) z &\text{if } \text{ } & a \le r \le b\\ &\bigg(-H_0 + \frac{\alpha_1}{r^3}\bigg) z &\text{if } \text{ } & r \ge b, \end{aligned} \right. \end{equation}

\begin{equation} \vec{\nabla} \Psi = \left\{ \begin{aligned} &\delta_1 \hat{k} &\text{if } \text{ } & r \le a\\ &-3\gamma_1\frac{xz}{r^5}\hat{i} -3\gamma_1\frac{yz}{r^5}\hat{j} + \bigg(\beta_1 + \gamma_1\frac{1}{r^3} - 3\gamma_1\frac{z^2}{r^5} \bigg)\hat{k} &\text{if } \text{ } & a \le r \le b\\ &-3\alpha_1\frac{xz}{r^5}\hat{i} -3\alpha_1\frac{yz}{r^5}\hat{j} + \bigg(-H_0 + \alpha_1\frac{1}{r^3} - 3\alpha_1\frac{z^2}{r^5} \bigg)\hat{k} &\text{if } \text{ } & r \ge b, \end{aligned} \right. \end{equation}

\[ \vec{H} = - \vec{\nabla} \Psi, \]

\[ \vec{B} = - \mu \vec{\nabla} \Psi, \]

where

\[ r = \sqrt{x^2 + y^2 +z^2}, \]

\[ \mu_r = \frac{\mu_1}{\mu_0}, \]

\begin{equation} \begin{aligned} &\Omega = \frac{(\mu_r - 1)}{(\mu_r + 2)} \frac{a^3}{b^3}, \\ &\gamma_1 = \frac{-3b^3 H_0 \Omega}{(2\mu_r+1)-2 (\mu_r-1) \Omega}, \\ &\beta_1 = \frac{(2\mu_r+1) \gamma_1}{(\mu-1)a^3}, \\ &\alpha_1 = \frac{-b^3 H_0+2\mu_r\gamma_1-\mu_rb^3\beta_1}{2}, \\ &\delta_1 = \frac{\mu_r a^3\beta_1-2\mu_r\gamma_1}{a^3}. \end{aligned} \end{equation}

The mesh of the problem domain is shown in the figure below.

The structure of the mesh is similar to the two-dimensional mesh described above. The main elements of the mesh are listed below.

The sphere of a radius $a$. It delineates the inner space of the shield. The permeability of the material inside the sphere equals $\mu_0$.
The spherical shell with an inner and outer radii of $a$ and $b$, respectively. It represents the shield. The permeability in this region equals $\mu_1$.
The region between the sphere of a radius $b$ and the cube with a side of $2d_3$. This region represents the space outside the shield. The permeability in this region equals $\mu_0$.
The cube with a side of $2d_2$. This cube delineates the local mesh. The error norms are computed inside the local mesh. The error norms outside the local mesh are discarded.
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 of constructing a mesh that has spherical interfaces and contains no tetrahedral cells: the cells around the origin must be hexahedral.
The cube with a side of $2d_3$ represents the infinity.

The H-field induced by the magnetization of the magnetic material of the shield vanish at infinity. Therefore, the total H-field at infinity equals the applied H-field. Thus, the Dirichlet boundary condition

\[ \Psi = -H_0 z \]

where $H_0 = \dfrac{1}{\mu_0} B_0$, can be applied on all facets of the boundary. This corresponds to the "DirichletOnly" mode of the program. On the other hand, the applied H-field runs parallel the $z$ axis. Therefore, its normal component on the four vertical facets of the boundary equals zero. Consequently, we can apply the Dirichlet boundary condition on the facets at $z=-d_3$ and $z=d_3$, and the homogeneous Neumann boundary condition,

\[ \mu \hat{n}\cdot\vec{\nabla} \Psi = 0, \]

on the four vertical facets. The program in the mode "DirichletNeumann" does exactly this. The problem domain shown above is truncated. That is, the boundary of the problem domain is at a finite distance from the origin. The domain of the initial problem is unbounded. For this reason, application of the Dirichlet boundary condition or a combination of the Dirichlet and Neumann boundary conditions will lead to a truncation error. We would like to observe the convergence rate without the truncation error. In general, the truncation error can be driven to be insignificantly small by increasing the parameter $d_3$, see the figure above. This, however, will increase the simulation time. Instead, we will keep the size of the problem domain relatively small and use the exact boundary conditions. That is, we will apply the Dirichlet boundary condition to all facets of the boundary and use the known exact solution (instead of $ -H_0 z $) to compute the right-hand side $\eta$. This will allow us to switch off the truncation error and observe the convergence rates free of the truncation error. This is what program does if switched into the "Exact" mode. The user can switch the program between the three modes by setting Setting::type_of_bc.

The boundary value problem

The following boundary value problem is solved in the "DirichletOnly" and "Exact" modes:

\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \mu \vec{\nabla} \Psi \big)=0 & \text{in} & \Omega & \text{(i)}\\ \text{(e)} &\Psi = \eta & \text{on} & \Gamma_{D1} & \text{(ii)}\\ \text{(e)} &\Psi_{+} = \Psi_{-} & \text{on} & \Gamma_{I1} & \text{(iii)}\\ \text{(n)}&\mu_{+}\hat{n}\cdot\vec{\nabla}\Psi_{+}-\mu_{-}\hat{n}\cdot\vec{\nabla}\Psi_{-}=0& \text{on} & \Gamma_{I1}&\text{(iv)}\\ \text{(e)} &\Psi_{+} = \Psi_{-} & \text{on} & \Gamma_{I2} & \text{(v)}\\ \text{(n)}&\mu_{+}\hat{n}\cdot\vec{\nabla}\Psi_{+}-\mu_{-}\hat{n}\cdot\vec{\nabla}\Psi_{-}=0& \text{on} & \Gamma_{I2}&\text{(vi)}. \end{array} \end{equation}

The following boundary value problem is solved in the "DirichletNeumann" mode:

\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \mu \vec{\nabla} \Psi \big)=0 & \text{in} & \Omega & \text{(i)}\\ \text{(e)} &\Psi = \eta & \text{on} & \Gamma_{D1} & \text{(ii)}\\ \text{(n)}&\mu\hat{n}\cdot\vec{\nabla}\Psi = 0 &\text{on}&\Gamma_{R1}&\text{(iii)}\\ \text{(e)} &\Psi_{+} = \Psi_{-} & \text{on} & \Gamma_{I1} & \text{(iv)}\\ \text{(n)}&\mu_{+}\hat{n}\cdot\vec{\nabla}\Psi_{+}-\mu_{-}\hat{n}\cdot\vec{\nabla}\Psi_{-}=0& \text{on} & \Gamma_{I1}&\text{(v)}\\ \text{(e)} &\Psi_{+} = \Psi_{-} & \text{on} & \Gamma_{I2} & \text{(vi)}\\ \text{(n)}&\mu_{+}\hat{n}\cdot\vec{\nabla}\Psi_{+}-\mu_{-}\hat{n}\cdot\vec{\nabla}\Psi_{-}=0& \text{on} & \Gamma_{I2}&\text{(vii)}. \end{array} \end{equation}

In "DirichletOnly" and "DirichletNeumann" modes the parameter $\eta$ is evaluated as

\[ \eta = -H_0 x \]

in the two-dimensional version of the experiment and as

\[ \eta = -H_0 z, \]

in the three-dimensional version of the experiment. In the "Exact" mode the parameter $\eta$ is computed by invoking the exact analytical solutions listed above.

First, a class template (SolverSLDI) derived from StaticScalarSolver::Solver is used to obtain the numerical solution, $\Psi$. After that, the numerical solution is fed to objects derived from the class templates StaticScalarSolver::ProjectPSItoH and StaticScalarSolver::ProjectPSItoB to calculate the auxiliary vector field $\vec{H}$ and magnetic field $\vec{B}$ as

\[ \vec{H} = - \vec{\nabla} \Psi \]

and

\[ \vec{B} = - \mu \vec{\nabla} \Psi, \]

respectively.

The program

This numerical experiment is implemented in accordance with the base structure described in here. The build process generates two executable files: sld-i-square and sld-i-cube. That is, one executable file for the two-dimensional version of the experiment, sld-i-square, and one executable file for the three-dimensional version of the experiment, sld-i-cube. To rebuild them, change into sld-i/build/Release directory and execute the following:

user@computer .../sld-i/build/Release$ ./clean

user@computer .../sld-i/build/Release$ ./build

Then the experiment needs to be executed again. To do so, change into the sld-i/bin/Release directory and execute run-all script there,

user@computer .../sld-i/build/Release$ cd ../../bin/Release

user@computer .../sld-i/bin/Release$ ./run-all

This will generate vtu files in the ssol-iii/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 sld-i/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into sld-i/gmsh directory and executing the following:

user@computer .../sld-i/gmsh$ ./clean

user@computer .../sld-i/gmsh$ ./build

This will generate a set of globally refined meshes in sld-i/gmsh/data.

The SettingsSLDI 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 program can be switched between three modes of operation, DirichletOnly, DirichletNeumann, and Exact, by assigning a proper value to Settings::type_of_bc.

Simulation results

Both, two- and three- dimensional versions of the experiment were executed under the following conditions: $a = 0.2[m]$, $b=0.4[m]$, $d_2=0.8[m]$, $d_3=2.0[m]$, $\mu_1 = 4.0 \mu_0$. All convergence tables were generated in the "Exact" mode, see above.

Two-dimensional experiment

Total magnetic scalar potential

The figure below illustrates a contour plot of the total magnetic scalar potential, $\Psi$, the output of an object of the SolverSLDI type.

The corresponding convergence table is presented below.

p	r	cells	dofs	$\\|e\\|_{L^2}$	$\alpha_{L^2}$	$\\|e\\|_{H^1}$	$\alpha_{H^1}$
1	10	2916	2953	2.55e-04	-	1.98e-02	-
1	11	3600	3641	2.07e-04	1.99	1.78e-02	0.97
1	12	4356	4401	1.71e-04	2.00	1.63e-02	0.97
1	13	5184	5233	1.43e-04	2.01	1.49e-02	0.98
2	10	2916	11737	8.04e-06	-	8.41e-04	-
2	11	3600	14481	5.89e-06	2.95	6.84e-04	1.96
2	12	4356	17513	4.44e-06	2.96	5.67e-04	1.96
2	13	5184	20833	3.43e-06	2.96	4.78e-04	1.97
3	10	2916	26353	2.63e-07	-	3.74e-05	-
3	11	3600	32521	1.75e-07	3.88	2.75e-05	2.93
3	12	4356	39337	1.22e-07	3.79	2.08e-05	2.94
3	13	5184	46801	8.89e-08	3.63	1.61e-05	2.95

The following notations were used in the table header:

p - the degree of the interpolating Lagrange polynomials that constitute the shape functions.
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.

Auxiliary H-field

The auxiliary H-field was calculated by feeding the numerically calculated total magnetic scalar potential, $\Psi$, to an object of the type StaticScalarSolver::ProjectPSItoH. The figure below illustrates the calculated auxiliary H-field.

The tangential component of the H-field is continuous on the interfaces between dissimilar materials. The mesh is made such that all interfaces between dissimilar materials are delineated by the faces of the mesh cells. That is, no interface cuts through a cell. Consequently, we can use the Nedelec finite elements to model the H-field as they provide the desired behavior of the approximated vector field: the tangential component is guaranteed to be continuous on the faces of all cells (including the faces that constitute the interfaces) and the normal component can be discontinuous if necessary. In other words, the choice of the Nedelec finite elements enforces the continuity of the tangential component of the auxiliary H-field on the interfaces between dissimilar materials.

The corresponding convergence table is presented below.

p	r	cells	dofs	$\\|e\\|_{L^2}$	$\alpha_{L^2}$
1	10	2916	5868	1.98e-02	-
1	11	3600	7240	1.78e-02	0.97
1	12	4356	8756	1.63e-02	0.97
1	13	5184	10416	1.49e-02	0.98
2	10	2916	23400	8.41e-04	-
2	11	3600	28880	6.84e-04	1.96
2	12	4356	34936	5.67e-04	1.96
2	13	5184	41568	4.78e-04	1.97
3	10	2916	52596	3.74e-05	-
3	11	3600	64920	2.75e-05	2.93
3	12	4356	78540	2.08e-05	2.94
3	13	5184	93456	1.61e-05	2.95

The following notations were used in the table header:

p - the degree of the FE_Q finite elements that were used to model the total magnetic scalar potential, $\Psi$. The degree of the FE_Nedelec finite elements that were used to model $\vec{H}$ field equals $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.

Magnetic field

The magnetic field, $\vec{B}$, was calculated by feeding the numerically calculated total magnetic scalar potential, $\Psi$, to an object of the type StaticScalarSolver::ProjectPSItoB. The figure below illustrates the calculated magnetic field.

The normal component of the magnetic field is continuous on the interfaces between dissimilar materials. The mesh is made such that all interfaces between dissimilar materials are delineated by the faces of the mesh cells. That is, no interface cuts through a cell. Consequently, we can use the Raviart-Thomas finite elements to model magnetic field as they provide the desired behavior of the approximated vector field: the normal component is guaranteed to be continuous on the faces of all cells (including the faces that constitute the interfaces) and the tangential component can be discontinuous if necessary. In other words, the choice of the Raviart-Thomas finite elements enforces the continuity of normal component of the magnetic field on the interfaces between dissimilar materials.

The corresponding convergence table is presented below.

p	r	cells	dofs	$\\|e\\|_{L^2} / \mu_0$	$\alpha_{L^2}$
1	10	2916	5868	2.77e-02	-
1	11	3600	7240	2.49e-02	1.04
1	12	4356	8756	2.25e-02	1.03
1	13	5184	10416	2.06e-02	1.03
2	10	2916	23400	1.25e-03	-
2	11	3600	28880	1.02e-03	1.94
2	12	4356	34936	8.44e-04	1.94
2	13	5184	41568	7.12e-04	1.95
3	10	2916	52596	4.78e-05	-
3	11	3600	64920	3.49e-05	2.97
3	12	4356	78540	2.63e-05	2.98
3	13	5184	93456	2.03e-05	2.98

The following notations were used in the table header:

p - the degree of the FE_Q finite elements that were used to model the total magnetic scalar potential, $\Psi$. The degree of the FE_RaviartThomas finite elements that were used to model the magnetic field, $\vec{B}$, equals $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} / \mu_0$ - the normalized $L^2$ error norm.
$\|e\|_{H^1}$ - the $H^1$ error norm.
$\alpha_{L^2}$ - the order of convergence of the normalized $L^2$ error norm.
$\alpha_{H^1}$ - the order of convergence of the $H^1$ error norm.