 |
Logbook
(07-04-2025)
Static problems
|
|
|
Logbook
(07-04-2025)
Static problems
|
|
Classes:
1) BatchMMSVTI
2) SolverMMSVTI_T
3) SolverMMSVTI_A
4) ExactSolutionMMSVTI_Jf
5) ExactSolutionMMSVTI_B
6) DirichletBC_MMSVTI_T
7) DirichletBC_MMSVTI_A
8) SettingsMMSVTI
Files:
1) mms-vt-i/src/main.cpp
2) mms-vt-i/include/solver.hpp
3) mms-vt-i/src/solver.cpp
4) mms-vt-i/include/exact_solution.hpp
5) mms-vt-i/src/exact_solution.cpp
6) mms-vt-i/src/static_vector_input.cpp
7) mms-vt-i/include/settings.hpp
The purpose of this numerical experiment is to test the StaticVectorSolver::Solver2 class template in the three-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 three spatial dimensions (dim=3) these two solvers differ from each other only at three points:
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 three-dimensional version of the mms-v/ numerical experiment. The difference is in how the vector field on the right-hand side of the partial differential equation, \(\vec{T}\), is fed to the solver. In the mms-v/ numerical experiment \(\vec{T}\) is fed to the solver in a form of an analytical expression. In this numerical experiment \(\vec{T}\) is computed numerically in a preprocessing step and then fed to the solver that computes \(\vec{A}\).
Another difference between the mms-v/ numerical experiment and this numerical experiment is in the amount of spatial dimensions in the solved problems. In the mms-v/ numerical experiment both, two- and three-dimensional problems are considered. In this numerical experiment - only three-dimensional. Test of the StaticVectorSolver::Solver2 class template in the two-dimensional space is done in the mms-v-ii/ 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 expressions for the volume free-current density, \(\vec{J}_M\), and other 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:
The last step, i.e., converting \(\vec{A}\) into the magnetic field \(\vec{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 \(\vec{B}\) we compute an error as
\[ e = | \vec{B}_M - \vec{B} |, \]
where
\[ \vec{B}_M = \vec{\nabla} \times \vec{A}_M \]
is the manufactured magnetic field. Then we evaluate the \(L^2\) error norm as
\[ ||e||_{L^2} = \sqrt{\iiint_{\Omega} e^2 dV}. \]
The three simulation steps above, i.e., computing \(\vec{T}\), \(\vec{A}\), and \(\vec{B}\), are mandatory. We add to them one optional simulation for checking if the numerically calculated current vector potential, \(\vec{T}\), has been calculated correctly. As \(\vec{T}\) is not gauged explicitly, we can evaluate it in the same way as we evaluate \(\vec{A}\). That is, we convert the numerically calculated current vector potential back to the volume free-current density,
\[ \vec{J}_{f} = \vec{\nabla} \times \vec{T}, \]
compute an error
\[ e = | \vec{J}_{M} - \vec{J}_{f} |, \]
where
\[ \vec{J}_{M} = \vec{\nabla}\times\bigg(\frac{1}{\mu}\vec{\nabla}\times\vec{A}_M\bigg), \]
is the manufactured free-current density. Then we evaluate the \(L^2\) error norm as
\[ ||e||_{L^2} = \sqrt{\iiint_{\Omega} e^2 dV}. \]
The four simulation steps (stages) discussed above are illustrated by means of the Bossavit's diagrams in the following figure.
As implied by the Bossavit's diagrams above, see also this discussion, the two 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, the magnetic field, \(\vec{B}\), and the free-current density, \(\vec{J}_f\), belong to the \(H(\text{div})\) function space, and, thus, must be modeled by the FE_RaviartThomas finite elements.
We will repeat the calculations of the \(L^2\) error norms on a family of progressively refined meshes. The \(L^2\) error norms are 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 as
\[ \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}. \]
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}\times\bigg(\frac{1}{\mu}\vec{\nabla}\times\vec{A}_M\bigg) = \vec{\nabla}\times\bigg(\frac{1}{\mu}\vec{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} \]
where
\[ B_M = \cos{(kx)} + \cos{(ky)}. \]
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} \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), \]
where
\[ \vec{T}_M = \frac{1}{\mu} \vec{\nabla}\times\vec{A}_M \]
is the manufactured current vector potential.
Two problem domains are used. The first problem domain is shaped as a cube. The second problem domain is shaped as a sphere. 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.
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) = \vec{\nabla}\times\vec{J}_M & \text{in} & \Omega & \text{(i)},\\ \text{(e)}& \hat{n}\times\vec{T} = \hat{n}\times\vec{T}_M & \text{on} & \Gamma_{D1} & \text{(ii)}, \\ \text{(n)}& \hat{n}\times\bigg(\vec{\nabla}\times\vec{T}\bigg) = \hat{n} \times \vec{J}_M & \text{on} & \Gamma_{R1} & \text{(iii)}, \\ \end{array} \end{equation}
where \(\vec{T}_M\) is the manufactured current vector potential and \(\vec{J}_M\) is the manufactured free-current density.
\[ \vec{J}_{f} = \vec{\nabla} \times \vec{T}. \]
Then \(\vec{J}_{f}\) is compared to the manufactured free-current density, \(\vec{J}_M\), by computing the \(L^2\) error norm.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}\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}
where \(\vec{T}\) is the current vector potential computed at the 0-th stage.
\[ \vec{B} = \vec{\nabla} \times \vec{A}. \]
Then \(\vec{B}\) is compared to the manufactured magnetic field, \(\vec{B}_M\), by computing the \(L^2\) error norm.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 were used:
The build process generates two executable files: mms-vt-i-cube, mms-vt-i-sphere. To rebuild them change into mms-vt-i/build/Release directory and execute the following:
Then all executable files must be executed again. This can be done by changing into the mms-vt-i/bin/Release directory and executing run-all script there,
This will generate vtu files in the mms-vt-i/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-i/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into mms-vt-i/gmsh directory and executing the following:
This will generate a set of refined meshes in mms-vt-i/gmsh/data.
The SettingsMMSVTI 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.
The simulations described below have been done with a setting \(k=\pi\). Recall, we do not gauge vector potentials explicitly. For this reason, presenting simulation results in terms of \(\vec{T}\) and \(\vec{A}\) does not make any sense: the solutions being expressed in terms of vector potentials are not unique. Instead, we present the results in terms of the curls of the vector potentials, i.e., in terms of \(\vec{J}_{f}\) and \(\vec{B}\).
The figure below illustrates the free-curent density calculated at the 1st stage, \(\vec{J}_{f}\), the output of an object of the StaticVectorSolver::ProjectTtoJ type. The vector field \(\vec{J}_f\) in the middle plot is perpendicular to the plane of cross section and, thus, invisible.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 0 | 5 | 5832 | 19494 | 5.73e+05 | - |
| 0 | 6 | 10648 | 34914 | 4.69e+05 | 1.00 |
| 0 | 7 | 17576 | 56862 | 3.97e+05 | 1.00 |
| 0 | 8 | 27000 | 86490 | 3.44e+05 | 1.00 |
| 1 | 5 | 5832 | 147852 | 3.03e+04 | - |
| 1 | 6 | 10648 | 267300 | 2.03e+04 | 1.99 |
| 1 | 7 | 17576 | 438204 | 1.45e+04 | 2.00 |
| 1 | 8 | 27000 | 669780 | 1.09e+04 | 2.00 |
| 2 | 5 | 5832 | 490050 | 1.26e+03 | - |
| 2 | 6 | 10648 | 888822 | 6.92e+02 | 2.99 |
| 2 | 7 | 17576 | 1460394 | 4.19e+02 | 3.00 |
| 2 | 8 | 27000 | 2235870 | 2.73e+02 | 3.00 |
The following notations were used in the table header:
The figure below illustrates the magnetic field, \(\vec{B}\), calculated at the 3rd stage, the output of an object of the StaticVectorSolver::ProjectAtoB type. The vector field \(\vec{B}\) in the right plot is perpendicular to the plane of cross section and, thus, invisible.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 0 | 5 | 5832 | 19494 | 2.86e-01 | - |
| 0 | 6 | 10648 | 34914 | 2.34e-01 | 1.00 |
| 0 | 7 | 17576 | 56862 | 1.98e-01 | 1.00 |
| 0 | 8 | 27000 | 86490 | 1.71e-01 | 1.00 |
| 1 | 5 | 5832 | 147852 | 1.29e-02 | - |
| 1 | 6 | 10648 | 267300 | 8.60e-03 | 2.00 |
| 1 | 7 | 17576 | 438204 | 6.16e-03 | 2.00 |
| 1 | 8 | 27000 | 669780 | 4.62e-03 | 2.00 |
| 2 | 5 | 5832 | 490050 | 3.79e-04 | - |
| 2 | 6 | 10648 | 888822 | 2.08e-04 | 3.00 |
| 2 | 7 | 17576 | 1460394 | 1.26e-04 | 3.00 |
| 2 | 8 | 27000 | 2235870 | 8.19e-05 | 3.00 |
The following notations were used in the table header:
The figure below illustrates the free-current density calculated at the 1st stage, \(\vec{J}_f\), the output of an object of the StaticVectorSolver::ProjectTtoJ type. The vector field \(\vec{J}_f\) in the middle plot is perpendicular to the plane of cross section and, thus, invisible.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 0 | 5 | 3584 | 11176 | 5.46e+05 | - |
| 0 | 6 | 7000 | 21650 | 4.39e+05 | 0.98 |
| 0 | 7 | 12096 | 37212 | 3.66e+05 | 0.99 |
| 0 | 8 | 19208 | 58870 | 3.14e+05 | 0.99 |
| 1 | 5 | 3584 | 87632 | 5.20e+04 | - |
| 1 | 6 | 7000 | 170500 | 3.45e+04 | 1.84 |
| 1 | 7 | 12096 | 293880 | 2.47e+04 | 1.83 |
| 1 | 8 | 19208 | 465836 | 1.86e+04 | 1.82 |
| 2 | 5 | 3584 | 293880 | 5.02e+03 | - |
| 2 | 6 | 7000 | 572550 | 2.79e+03 | 2.63 |
| 2 | 7 | 12096 | 987732 | 1.73e+03 | 2.63 |
| 2 | 8 | 19208 | 1566642 | 1.15e+03 | 2.63 |
The following notations were used in the table header:
The figure below illustrates the magnetic field, \(\vec{B}\), calculated at the 3rd stage, the output of an object of the StaticVectorSolver::ProjectAtoB type. The vector field \(\vec{B}\) in the right plot is perpendicular to the plane of cross section and, thus, invisible.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 0 | 5 | 3584 | 11176 | 2.10e-01 | - |
| 0 | 6 | 7000 | 21650 | 1.67e-01 | 1.03 |
| 0 | 7 | 12096 | 37212 | 1.39e-01 | 1.02 |
| 0 | 8 | 19208 | 58870 | 1.19e-01 | 1.02 |
| 1 | 5 | 3584 | 87632 | 1.66e-02 | - |
| 1 | 6 | 7000 | 170500 | 1.09e-02 | 1.89 |
| 1 | 7 | 12096 | 293880 | 7.76e-03 | 1.87 |
| 1 | 8 | 19208 | 465836 | 5.83e-03 | 1.85 |
| 2 | 5 | 3584 | 293880 | 1.27e-03 | - |
| 2 | 6 | 7000 | 572550 | 7.09e-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