 |
Logbook
(07-04-2025)
Static problems
|
|
|
Logbook
(07-04-2025)
Static problems
|
|
Classes:
1) BatchCVPI
2) SolverCVPI
3) ExactSolutionCVPI_Jf
4) DirichletBC_CVPI
5) SettingsCVPI
Files:
1) cvp-i/src/main.cpp
2) cvp-i/include/solver.hpp
3) cvp-i/src/solver.cpp
4) cvp-i/include/exact_solution.hpp
5) cvp-i/src/exact_solution.cpp
6) cvp-i/src/static_vector_input.cpp
7) cvp-i/include/settings.hpp
In this numerical experiment we will convert the free-current density, \(\vec{J}_f\), shown below into the current vector potential, \(\vec{T}\).
The compatibility condition for the curl-curl equation,
\[ \vec{\nabla}\times\bigg(\frac{1}{\mu}\vec{\nabla}\times\vec{A}\bigg) = \vec{J}_f, \]
requires the free-current density, \(\vec{J}_f\), to be derived as a curl of a current vector potential,
\[ \vec{J}_f = \vec{\nabla}\times\vec{T} \]
(see Section 5.5. Compatibility conditions). Then the initial curl-curl equation becomes
\[ \vec{\nabla}\times\bigg(\frac{1}{\mu}\vec{\nabla}\times\vec{A}\bigg) = \vec{\nabla}\times\vec{T}. \]
In order to solve the last equation, one needs to convert the free-current density, \(\vec{J}_f\), given by the formulation of the problem into a current vector potential, \(\vec{T}\). The most straightforward way to do so is to solve the following partial differential equation:
\[ \vec{\nabla}\times\bigg(\vec{\nabla}\times\vec{T}\bigg) = \vec{\nabla}\times\vec{J}_f. \]
If the free-current density, \(\vec{J}_f\), is given in a form of a closed-form analytical expression, the last equation can be solved with a help of StaticVectorSolver::Solver1 class template. This numerical experiment illustrates how to do that.
First, we will solve the last partial differential equation for \(\vec{T}\) assuming that the free-current density, \(\vec{J}_f\), is given in a form of analytical expression. Second, we will convert the numerically calculated \(\vec{T}\) back into \(\vec{J}_{f}\) for evaluation. The Bossavit's diagrams below illustrate these two steps.
Third, we will compare the numerically calculated \(\vec{J}_{f}\) with the initial closed-form analytical expression by computing the \(L^2\) error norm. Note that we will not gauge the current vector potential explicitly. Its conservative part will remain unknown. For this reason, it does not make much sense to visualize it or evaluate it by computing error norms. For this reason, we will evaluate the computed current vector potential, \(\vec{T}\), indirectly by computing \(L^2\) error norm of the free-current density, \(\vec{J}_f\), derived from it. As per usual, we will repeat the simulations on a family of progressively refined meshes and observe the behavior of the \(L^2\) error norm.
As implied by the Bossavit's diagrams above, see also this discussion, the current vector potential, \(\vec{T}\), belongs to the \(H(\text{curl})\) function space, and, thus, must be modeled by the FE_Nedelec finite elements. Similarly, the free-current density, \(\vec{J}_f\), belongs to the \(H(\text{div})\) function space, and, thus, must be modeled by the FE_RaviartThomas finite elements.
Note that this numerical experiment is carried out in three spatial dimensions. The current vector potential calculated in this numerical experiment, \(\vec{T}\), can be used in three-dimensional problems only. To solve for the current potential that can be used in two-dimensional problems, \(T\), (it is always an out-of-plane vector, i.e., a scalar) one needs to invoke the div-grad equation. How to do this is illustrated in the Current vector potential (cvp-ii) numerical experiment.
The problem domain is assumed to be a ball concentric with the origin. The radius of the ball equals \( d_2 \). The boundary of the ball represents infinity. We apply the homogeneous Dirichlet boundary condition,
\[ \hat{n} \times \vec{T} = 0. \]
The volume free-current density, \( \vec{J}_f \), is assumed to be confined to a spherical shell concentric with the origin. The inner and outer radii of the shell equal \( a \) and \( b \), respectively. The figure below illustrates the problem domain.
The free-current density can be expressed as
\[ \vec{J}_f = \left\{ \begin{aligned} & 0 &&\text{if }&& r < a\\ & -y\hat{i} + x\hat{j} && \text{if } && a < r < b\\ & 0 &&\text{if }&& r > b, \end{aligned} \right. \]
where
\[ r = \sqrt{x^2+y^2+z^2}. \]
Recall that we do not explicitly gauge the current vector potential, \(\vec{T}\). For this reason, we do not specify a closed-form analytical equation for it as its conservative part is unknown.
The following boundary value problem is solved in this numerical experiment:
\begin{equation} \begin{array}{lrcll} \text{ } & \vec{\nabla}\times\bigg(\vec{\nabla}\times\vec{T}\bigg) = \vec{\nabla}\times\vec{J}_f & \text{in} & \Omega & \text{(i)},\\ \text{(e)}&\hat{n} \times \vec{T} = 0 & \text{on} & \Gamma_{D1} & \text{(ii)}. \\ \end{array} \end{equation}
To verify the solution, the numerically computed current vector potential, \(\vec{T}\), is converted back to the free-current density as
\[ \vec{J}_{f} = \vec{\nabla} \times \vec{T}, \]
and compared to the free-current density, \(\vec{J}_f\), given by the closed-form analytical expression above.
This experiment is implemented in accordance with the base code structure. The build process generates one executable file: cvp-i-linear. To rebuild it, change into cvp-i/build/Release directory and execute the following:.
Then the experiment needs to be executed again. To do so, change into the cvp-i/bin/Release directory and run the executable file,
This will generate vtu files in the cvp-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.
Note that executable files require a set of meshes to be present in the cvp-i/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into cvp-i/gmsh directory and executing the following:
This will generate a set of globally refined meshes in cvp-i/gmsh/data.
The SettingsCVPI 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.
It worth mentioning that the free-current density, \(\vec{J}_f\), equals zero on the boundary by the definition of the problem. For this reason we can neglect the boundary integral \(I_{b3-2}\) in the Recipe for static vector solver in 3D (current vect. potential). Consequently, we pass the input parameter type_of_pde_pde = 2 to the constructor of the StaticVectorSolver::Solver1:
The experiment was executed with the following settings: \(a = 0.3[m]\), \(b = 0.5[m]\), and \(d_2=1.0[m]\).
Two two-dimensional slices of the numerically calculated free-current density are presented in the figure below. The vectors of the \(\vec{J}_{f}\) field in the left plot are perpendicular to the slice and, thus, invisible.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 0 | 12 | 25289 | 76230 | 1.03e-02 | - |
| 0 | 13 | 32832 | 98928 | 9.47e-03 | 1.00 |
| 0 | 14 | 41743 | 125736 | 8.74e-03 | 1.00 |
| 0 | 15 | 52136 | 156996 | 8.12e-03 | 1.00 |
| 1 | 12 | 25289 | 608388 | 3.29e-05 | - |
| 1 | 13 | 32832 | 789696 | 2.75e-05 | 2.05 |
| 1 | 14 | 41743 | 1003860 | 2.33e-05 | 2.05 |
| 1 | 15 | 52136 | 1253616 | 2.01e-05 | 2.04 |
| 2 | 12 | 25289 | 2051676 | 7.94e-06 | - |
| 2 | 13 | 32832 | 2663280 | 6.12e-06 | 2.99 |
| 2 | 14 | 41743 | 3385746 | 4.82e-06 | 2.99 |
| 2 | 15 | 52136 | 4228308 | 3.87e-06 | 2.98 |
The following notations were used in the table header:
1.9.1