 |
Logbook
(07-04-2025)
Static problems
|
|
|
Logbook
(07-04-2025)
Static problems
|
|
Classes:
1) BatchMMSV
2) SolverMMSV
3) ExactSolutionMMSV_A
4) ExactSolutionMMSV_B
5) SettingsMMSV
Files:
1) mms-v/src/main.cpp
2) mms-v/include/solver.hpp
3) mms-v/src/solver.cpp
4) mms-v/include/exact_solution.hpp
5) mms-v/src/exact_solution.cpp
6) mms-v/src/static_vector_input.cpp
7) mms-v/include/settings.hpp
The purpose of this numerical experiment is to verify the code of the following three class templates
by means of the method of manufactured solutions. The class template StaticVectorSolver::Solver1 implements a solver for the static vector boundary value problem. The other two templates, StaticVectorSolver::ProjectAtoB and StaticVectorSolver::ProjectAxyToBz, convert the magnetic vector potential, \(\vec{A}\), into magnetic field \(\vec{B}\). This numerical experiment is similar to the Method of manufactured solutions, scalar potential (mms/) numerical experiment. There are, however, two significant differences between this experiment and Method of manufactured solutions, scalar potential (mms/).
The first difference is due to the fact that we will not gauge the magnetic vector potential explicitly. This means that the solution to the curl-curl equation will not be unique. It does not make any sense to speak about visualization or evaluation of computed magnetic vector potential, \(\vec{A}\), as its conservative part will remain unknown, and, as a consequence, we cannot write down the corresponding exact closed-form analytical expression. Instead, we will evaluate the quality of the calculated magnetic vector potential indirectly by converting it into magnetic field,
\[ \vec{B}= \vec{\nabla} \times \vec{A}, \]
and then calculating the \(L^2\) error norms by comparing \(\vec{B}\) with its exact closed-form analytical expression. The magnetic field will be unique even if the magnetic vector potential is not gauged. Consequently, it is possible to manufacture closed-form expression for the magnetic field. In such disposition, however, the numerically calculated magnetic field is a result of a collective effort of two algorithms:
We cannot dissociate errors made by the two algorithms. Consequently, we have to evaluate the total error made by both algorithms.
The second significant difference between this numerical experiment and Method of manufactured solutions, scalar potential (mms/) originates in the nature of the curl-curl equation. The cur-curl equation in three dimensions reads
\begin{equation} \vec{\nabla}\times\bigg(\frac{1}{\mu}\vec{\nabla}\times\vec{A}\bigg) = \vec{J}_f. \end{equation}
The curl-curl equation in two dimensions reads
\begin{equation} \vec{\nabla} \overset{V}{\times}\bigg(\frac{1}{\mu} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg) = \vec{J}_f. \end{equation}
The compatibility condition for the curl-curl equation requires that the free-current density, \(\vec{J}_f\), is derived as a curl of a current vector potential, \(\vec{T}\). Then the curl-curl equation in three and two dimensions can be expressed as
\begin{equation} \vec{\nabla}\times\bigg(\frac{1}{\mu}\vec{\nabla}\times\vec{A}\bigg)= \vec{\nabla}\times\vec{T} \end{equation}
and
\begin{equation} \vec{\nabla} \overset{V}{\times}\bigg(\frac{1}{\mu}\vec{\nabla}\overset{S}{\times}\vec{A}\bigg) =\vec{\nabla} \overset{V}{\times} T, \end{equation}
respectively. Strictly speaking, we cannot just place the sources, \(\vec{J}_f\), on the right-hand side of the curl-curl equation. On the contrary, placing the sources, \(\rho_f\), on the right-hand side of the div-grad equation in the Method of manufactured solutions, scalar potential (mms/) numerical experiment is acceptable.
This numerical experiment consists of the following. First, we write down an arbitrary closed-form expression and assume that this expression is the solution to the boundary value problem. Let us call this expression a manufactured solution and denote it by \(\vec{A}_M\). Second, we differentiate the manufactured solution several times to get closed-form expressions of the right-hand sides in the boundary value problem, i.e., \(\vec{T}\), \(\vec{G}\), and \(\vec{Q}\), as well as the expression for the manufactured magnetic field, \(\vec{B}_M\). We choose closed-form expressions for \(\mu\) and \(\gamma\) somewhat arbitrary. Third, we obtain a numerical solution, \(\vec{A}\), by feeding the closed-form expressions for \(\vec{T}\), \(\vec{G}\), \(\vec{Q}\), \(\mu\), and \(\gamma\) to the StaticVectorSolver::Solver1. Fourth, we convert \(\vec{A}\) to magnetic field, \(\vec{B}\), with a help of StaticVectorSolver::ProjectAtoB in 3D or StaticVectorSolver::ProjectAxyToBz in 2D. The last two steps, i.e., calculating \(\vec{A}\) and convetring it into \(\vec{B}\), in the three-dimensional version of the experiment can be expressed with a help of the Bossavit's diagrams as the following: 
The figure below expresses the same two calculations in the two-dimensional version of the experiment: 
Next, we calculate an error as a magnitude of a difference between the manufactured and numerical solutions,
\[ e = | \vec{B}_M-\vec{B} |. \]
Finally, we evaluate the \(L^2\) error norm as
\[ ||e||_{L^2} = \sqrt{\iiint_{\Omega} e^2 dV} \\ \]
and as
\[ ||e||_{L^2} = \sqrt{\iint_{\Omega} e^2 dS} \\ \]
in three- and two- dimensional versions of the experiment, respectively.
As implied by the Bossavit's diagrams above, see also this discussion, the magnetic vector potential, \(\vec{A}\), in three- and two- dimensional versions of the experiment belongs to the \(H(\text{curl})\) function space, and, thus, must be modeled by the FE_Nedelec finite elements. Similarly, the magnetic field, \(\vec{B}\), belongs to the \(H(\text{div})\) function space in the three- dimensional version of the experiment, and, thus, must be modeled by the FE_RaviartThomas finite elements. In the two-dimensional version of the experiment, the magnetic field belongs to the \(L_2\) function space, and, thus, must be modeled by the FE_DGQ finite elements.
We will repeat the calculations of the \(L^2\) error norm on a family of progressively refined meshes. The \(L^2\) error norm is supposed to converge to zero as the size of the mesh cell decreases. If it will converge to zero at any rate, we will conclude that the amount of bugs in the code is acceptable.
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,
\[ 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)}. \]
This equation will be used to calculate the \(L^2\) error norm. 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 current vector potential as
\[ T = \frac{1}{\mu} \vec{\nabla} \overset{S}{\times} \vec{A}_M = \frac{1}{\mu} B_M = \frac{\cos{(kx)} + \cos{(ky)} }{\mu_0 (x^2 + y^2 + 1)}. \]
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 + \gamma \hat{n} \overset{V}{\times} \bigg(\hat{n}\overset{S}{\times}\vec{A}_M\bigg). \]
This numerical experiment reuses the two-dimensional problem domains, i.e., the square and the circle, of the Method of manufactured solutions, scalar potential (mms/) numerical experiment. The meshes in this numerical experiments and in Method of manufactured solutions, scalar potential (mms/) numerical experiment are identical:
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,
\[ \vec{B}_M = \vec{\nabla} \times \vec{A}_M = \begin{vmatrix} \hat{i} & \hat{j} & \hat{k} \\ \dfrac{\partial}{\partial x}&\dfrac{\partial}{\partial y}&\dfrac{\partial}{\partial z} \\ -\dfrac{1}{k}\sin{(ky)} & \dfrac{1}{k}\sin{(kx)} & 0 \end{vmatrix} = \big( \cos{(kx)} + \cos{(ky)} \big) \hat{k}. \]
This equation will be used to calculate the \(L^2\) error norm. 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 current vector potential as
\[ \vec{T} = \frac{1}{\mu} \vec{\nabla} \times \vec{A}_M = \frac{1}{\mu} \vec{B}_M = \frac{\big( \cos{(kx)} + \cos{(ky)} \big) \hat{k} }{\mu_0 (x^2 + y^2 + 1)}. \]
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} \times \bigg(\vec{\nabla}\times\vec{A}_M\bigg) + \gamma \hat{n} \times \bigg(\hat{n}\times\vec{A}_M\bigg) = \hat{n} \times \vec{T}_M + \gamma \hat{n} \times \bigg(\hat{n}\times\vec{A}_M\bigg). \]
This numerical experiment reuses the three-dimensional problem domains, i.e., the cube and the sphere, of the Method of manufactured solutions, scalar potential (mms/) numerical experiment. The meshes in this numerical experiments and in Method of manufactured solutions, scalar potential (mms/) numerical experiment are identical:
The two-dimensional version of this experiment solves 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}
The three-dimensional version of this experiment solves the following boundary value problem:
\begin{equation} \begin{array}{lrcll} \text{ } & \vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A}\bigg) = \vec{\nabla}\times\vec{T} & \text{in} & \Omega & \text{(i)},\\ \text{(e)}&\hat{n} \times \vec{A} = \hat{n} \times \vec{G} & \text{on} & \Gamma_{D1} & \text{(ii)}, \\ \text{(n)}& \dfrac{1}{\mu}\hat{n}\times\bigg(\vec{\nabla}\times\vec{A}\bigg) + \gamma \hat{n} \times \bigg( \hat{n}\times\vec{A} \bigg) = \vec{Q} & \text{on} & \Gamma_{R1} & \text{(iii)}. \\ \end{array} \end{equation}
These boundary value problems are special cases of the general two- and three-dimensional static vector boundary value problems
This numerical experiment is implemented in accordance with the base structure described in here. The class template StaticVectorSolver::ProjectAtoB does not use any external coefficients. Therefore, the class templates StaticVectorSolver::Solver1 and StaticVectorSolver::ProjectAtoB can share the same set of class templates declared in static_vector_input.hpp. Consequently, the class templates StaticVectorSolver::Solver1 and StaticVectorSolver::ProjectAtoB in the three-dimensional version of the experiment are instantiated at the same stage, the default first stage. The same is true for the class templates StaticVectorSolver::Solver1 and StaticVectorSolver::ProjectAxyToBz in the two-dimensional version of the experiment.
The build process generates four executable files: mms-v-square, mms-v-circle, mms-v-cube, and mms-v-sphere. To rebuild them change into mms-v/build/Release directory and execute the following:
Then all executable files must be executed again. This can be done by changing into the mms-v/bin/Release directory and executing run-all script there,
This will generate various files in the mms-v/bin/Release/Data directory. Among the generated files there are vtu files. 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 on the circular and spherical meshes.
That executable files require a set of meshes to be present in the mms-v/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:
This will generate a set of refined meshes in mms-v/gmsh/data.
The SettingsMMSV 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.
All simulations have been done with a setting \(k=\pi\).
The figure below illustrates a plot of the calculated magnetic field, \(B\), the output of an object of the StaticVectorSolver::ProjectAxyToBz type.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 0 | 5 | 64 | 144 | 4.49e-01 | - |
| 0 | 6 | 100 | 220 | 3.60e-01 | 0.98 |
| 0 | 7 | 144 | 312 | 3.01e-01 | 0.99 |
| 0 | 8 | 196 | 420 | 2.58e-01 | 0.99 |
| 1 | 5 | 64 | 544 | 4.56e-02 | - |
| 1 | 6 | 100 | 840 | 2.93e-02 | 1.99 |
| 1 | 7 | 144 | 1200 | 2.04e-02 | 1.99 |
| 1 | 8 | 196 | 1624 | 1.50e-02 | 1.99 |
| 2 | 5 | 64 | 1200 | 3.03e-03 | - |
| 2 | 6 | 100 | 1860 | 1.56e-03 | 2.99 |
| 2 | 7 | 144 | 2664 | 9.02e-04 | 2.99 |
| 2 | 8 | 196 | 3612 | 5.68e-04 | 3.00 |
The following notations were used in the table header:
The figure below illustrates a plot of the calculated magnetic field, \(B\), the output of an object of the StaticVectorSolver::ProjectAxyToBz type.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 0 | 5 | 192 | 400 | 2.89e-01 | - |
| 0 | 6 | 300 | 620 | 2.31e-01 | 1.00 |
| 0 | 7 | 432 | 888 | 1.93e-01 | 1.00 |
| 0 | 8 | 588 | 1204 | 1.65e-01 | 1.00 |
| 1 | 5 | 192 | 1568 | 1.49e-02 | - |
| 1 | 6 | 300 | 2440 | 9.49e-03 | 2.03 |
| 1 | 7 | 432 | 3504 | 6.56e-03 | 2.02 |
| 1 | 8 | 588 | 4760 | 4.81e-03 | 2.02 |
| 2 | 5 | 192 | 3504 | 8.01e-04 | - |
| 2 | 6 | 300 | 5460 | 4.05e-04 | 3.06 |
| 2 | 7 | 432 | 7848 | 2.32e-04 | 3.05 |
| 2 | 8 | 588 | 10668 | 1.45e-04 | 3.04 |
The following notations were used in the table header:
The figure below illustrates a plot of the calculated magnetic field, \(B\), the output of an object of the StaticVectorSolver::ProjectAtoB type.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 0 | 5 | 512 | 1944 | 6.35e-01 | - |
| 0 | 6 | 1000 | 3630 | 5.10e-01 | 0.99 |
| 0 | 7 | 1728 | 6084 | 4.26e-01 | 0.99 |
| 0 | 8 | 2744 | 9450 | 3.65e-01 | 0.99 |
| 1 | 5 | 512 | 13872 | 6.45e-02 | - |
| 1 | 6 | 1000 | 26460 | 4.14e-02 | 1.99 |
| 1 | 7 | 1728 | 45000 | 2.88e-02 | 1.99 |
| 1 | 8 | 2744 | 70644 | 2.12e-02 | 1.99 |
| 2 | 5 | 512 | 45000 | 4.29e-03 | - |
| 2 | 6 | 1000 | 86490 | 2.20e-03 | 2.99 |
| 2 | 7 | 1728 | 147852 | 1.28e-03 | 2.99 |
| 2 | 8 | 2744 | 232974 | 8.04e-04 | 3.00 |
The following notations were used in the table header:
The figure below illustrates a plot of the calculated magnetic field, \(B\), the output of an object of the StaticVectorSolver::ProjectAtoB type.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 0 | 5 | 3584 | 11176 | 2.09e-01 | - |
| 0 | 6 | 7000 | 21650 | 1.67e-01 | 1.02 |
| 0 | 7 | 12096 | 37212 | 1.38e-01 | 1.01 |
| 0 | 8 | 19208 | 58870 | 1.18e-01 | 1.01 |
| 1 | 5 | 3584 | 87632 | 1.65e-02 | - |
| 1 | 6 | 7000 | 170500 | 1.09e-02 | 1.88 |
| 1 | 7 | 12096 | 293880 | 7.74e-03 | 1.86 |
| 1 | 8 | 19208 | 465836 | 5.82e-03 | 1.85 |
| 2 | 5 | 3584 | 293880 | 1.27e-03 | - |
| 2 | 6 | 7000 | 572550 | 7.07e-04 | 2.61 |
| 2 | 7 | 12096 | 987732 | 4.41e-04 | 2.60 |
| 2 | 8 | 19208 | 1566642 | 2.96e-04 | 2.59 |
The following notations were used in the table header:
1.9.1