 |
Logbook
(07-04-2025)
Static problems
|
|
|
Logbook
(07-04-2025)
Static problems
|
|
Classes:
1) BatchCBND
2) SolverCBND
3) ExactSolutionCBND_PHI
4) SettingsCBND
Files:
1) cbnd/src/main.cpp
2) cbnd/include/solver.hpp
3) cbnd/src/solver.cpp
4) cbnd/include/exact_solution.hpp
5) cbnd/src/exact_solution.cpp
6) cbnd/src/static_scalar_input.cpp
7) cbnd/include/settings.hpp
This numerical experiment allows observing the influence of a curved boundary of a problem domain on convergence rates. In the Method of manufactured solutions (mms/) numerical experiment high convergence orders were observed:
\[ \alpha_{L^2} \approx p+1 \]
and
\[ \alpha_{H^1} \approx p, \]
where \(p\) is the degree of the interpolating polynomials. These are the theoretical upper boundaries of the convergence orders. The convergence orders were so high due to the fact that we have used the manufactured solution to calculate the right-hand sides of the Dirichlet and Robin boundary conditions at exact locations of the support points. In a real-life simulation there is no manufactured solution. Consequently, we can calculate the right-hand sides of the boundary conditions at the exact locations of the support points only if the problem domain has a polygonal shape and mesh approximates the boundary of the problem domain exactly. The following figure illustrates this situation.
If the boundary of the problem domain is curved and the degree of the interpolating polynomials is \(p>1\), some of the support points are displaced from the boundary of the domain. As a result, the solver applies the boundary conditions at the support points that are displaced from the points at which the boundary conditions are specified. The following figure illustrates this situation.
The fact that a polygonal mesh cannot approximate a curved boundary of a problem domain exactly contributes to the total error of the simulation decreasing the convergence rate. It is well known (see [1], for instance) that if a boundary of a problem domain is shaped as a circle, the \(L^2\) convergence rate decreases to
\[ \alpha_{L^2} = 2, \]
which is lower than the theoretical upper bound for \( p > 1\).
Observation of this degradation of the convergence rate is the first goal of this numerical experiment. The second goal is to test the code implemented by the class StaticScalarSolver::Solver by the method of exact solutions. The method of exact solutions is identical to the Method of manufactured solutions (mms/) with one exception: the method of exact solutions utilizes closed-form solutions to textbook problems instead of manufactured solutions.
The following textbook problem is solved in the two-dimensional version of this numerical experiment. (Section 6.1.1. Two coaxial cylindrical tubes)
The system consists of two infinitely long conducting coaxial cylindrical tubes. The space between the tubes is empty. The first tube is at a potential of \(\Phi=\Phi_0\) Volts. The second tube is grounded, so the electric field exists only in the space between the tubes.
We would like to calculate the electric potential.
As soon as the system exhibits translation symmetry, a two-dimensional planar problem domain can be used to describe the system. We remove the conductors from the problem domain and model them by setting up Dirichlet boundary conditions on the surfaces that delineate them. This problem can be easily solved by a number of analytical methods. The closed-form solutions are listed below.
\begin{equation} \Phi = \frac{\ln(b/r)}{\ln(b/a)} \Phi_0, \end{equation}
\begin{equation} \vec{\nabla} \Phi = -\frac{\Phi_0}{\ln(b/a)} \frac{1}{r}\hat{r}. \end{equation}
The two expressions above are used for estimating the \(L^2\) and \(H^1\) error norms.
The following textbook problem is solved in the three-dimensional version of this numerical experiment. (Section 6.1.2. Two concentric spherical shells)
The system consists of two conducting concentric spherical shells. The space between the shells is empty. The first shell is at a potential of \(\Phi = \Phi_0\) Volts. The second shell is grounded, so the electric field exists only in the space between the shells. The figure below illustrates the system and the problem domain.
We would like to calculate the electric potential.
Because the problem domain exhibits spherical symmetry, it can be reduced to a one-dimensional problem. However, we will use the three-dimensional problem domain as our intention is to test the solver in the three-dimensional space. We remove the conductors from the problem domain and model them by setting up Dirichlet boundary conditions on the surfaces that delineate them. This problem can be easily solved by a number of analytical methods. The closed-form solutions are listed below.
\begin{equation} \Phi = \frac{ab}{b - a} \Bigg(\frac{1}{r} - \frac{1}{b}\Bigg) \Phi_0, \end{equation}
\begin{equation} \vec{\nabla} \Phi = - \frac{ab}{b-a} \frac{\Phi_0}{r^2}\hat{r}. \end{equation}
The two expressions above are used for estimating the \(L^2\) and \(H^1\) error norms.
Both versions of this numerical experiment, the two- and three- dimensional, solve the following boundary value problem by the finite element method:
\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \epsilon_0 \vec{\nabla} \Phi \big)= 0 & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &\Phi = \eta_1 & \text{on} & \Gamma_{D1} & \text{(ii)},\\ \text{(e)} &\Phi = \eta_2 & \text{on} & \Gamma_{D2} & \text{(iii)}. \end{array} \end{equation}
This boundary value problem is a special case of the general static scalar boundary value problem. The exact values of the parameters \(\eta_1\) and \(\eta_2\) depend on the macro definition IS_BC_EXACT__ set in the CMakeLists.txt file. If IS_BC_EXACT__=0, the parameters take the following values
\[ \eta_1 = \Phi_0 \]
and
\[ \eta_2 = 0, \]
where \(\Phi_0\) is a constant. In the case IS_BC_EXACT__=1, both, \(\eta_1\) and \(\eta_2\), are calculated by evaluating the closed-form expressions for \(\Phi\) given above.
The CBND experiment is implemented in accordance with the base code structure. The build process generates four executable files: cbnd-ring, cbnd-shell, cbnd-ring-exact, and cbnd-shell-exact. That is, two executable files for the two- dimensional experiment: cbnd-ring, cbnd-ring-exact; and two executable files for the three- dimensional experiment: cbnd-shell, cbnd-shell-exact. To rebuild them change into cbnd/build/Release directory and execute the following:
Then all executable files must be executed again. This can be done by changing into the cbnd/bin/Release directory and executing run-all script there,
This will generate various files in the cbnd/bin/Release/Data directory. Among the generated files there are vtk files that can be viewed with a help of VisIt software package of the Lawrence Livermore National Laboratory. The data files in tex and txt format contain the convergence tables. The cbnd/bin/Release/Data directory also contains gpi files. They can be used to visualize the calculated potential along the \(x\) axis.
Note that executable files require a set of meshes to be present in the cbnd/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into cbnd/gmsh directory and executing the following:
This will generate a set of globally refined meshes in cbnd/gmsh/data.
The SettingsCBND 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 vtk files next to the numerical solution.
The suffix "exact" in the names of the executable files indicates that the Dirichlet boundary conditions on both boundaries were replaced with the exact solution. This replacement is done simply by replacing the functions in the stack of Dirichlet conditions as illustrated by the code snippet below.
As the result of this replacement, the values for the Dirichlet boundary conditions are calculated at exact locations. Consequently, the contribution of the boundary approximation error to the total error vanishes. Below, the Dirichlet boundary conditions calculated at exact locations will be called the exact boundary conditions. The exact boundary conditions allow switching off the influence of the curvature of the boundaries on the total simulation error.
The figure below illustrates the electric potential calculated for \(a=0.5 [m] \), \(b=1.0 [m] \), and \(\Phi_0 = 1 [V] \).
The table below illustrates convergence of the error norms with the Dirichlet boundary conditions applied on both boundaries.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) | \(\|e\|_{H^1}\) | \(\alpha_{H^1}\) |
|---|---|---|---|---|---|---|---|
| 1 | 15 | 784 | 840 | 2.14e-03 | - | 1.08e-01 | - |
| 1 | 16 | 900 | 960 | 1.86e-03 | 2.00 | 1.01e-01 | 1.00 |
| 1 | 17 | 1024 | 1088 | 1.64e-03 | 2.00 | 9.43e-02 | 1.00 |
| 1 | 18 | 1156 | 1224 | 1.45e-03 | 2.00 | 8.87e-02 | 1.00 |
| 2 | 15 | 784 | 3248 | 2.33e-03 | - | 2.06e-02 | - |
| 2 | 16 | 900 | 3720 | 2.03e-03 | 2.00 | 1.86e-02 | 1.50 |
| 2 | 17 | 1024 | 4224 | 1.79e-03 | 2.00 | 1.68e-02 | 1.50 |
| 2 | 18 | 1156 | 4760 | 1.58e-03 | 2.00 | 1.54e-02 | 1.50 |
| 3 | 15 | 784 | 7224 | 2.33e-03 | - | 2.00e-02 | - |
| 3 | 16 | 900 | 8280 | 2.03e-03 | 2.00 | 1.80e-02 | 1.50 |
| 3 | 17 | 1024 | 9408 | 1.79e-03 | 2.00 | 1.64e-02 | 1.50 |
| 3 | 18 | 1156 | 10608 | 1.58e-03 | 2.00 | 1.50e-02 | 1.50 |
The following notations were used in the table header:
The \(L^2\) error norm converges to zero with convergence order of
\[ \alpha_{L^2} = 2 \]
for all values of \(p\). This convergence order is lower than the convergence order we have observed in the Method of manufactured solutions (mms/) numerical experiment for \(p>1\). The \(H^1\) convergence order experiences a similar reduction. The reason for this reduction of the convergence orders is the curvature in the shape of the boundary of the problem domain. This assumption can be confirmed by the applying exact boundary condition on both boundaries. This can be done by running the cbnd-ring-exact executable file. The results are listed in the table below. The convergence orders in this case are the same as in the Method of manufactured solutions (mms/) numerical experiment. The last suggests that the curvature in the boundary of the problem domain is the reason for the reduction of the convergence orders in the case when the Dirichlet boundary conditions are applied on both boundaries in the usual way, i.e. \(\Phi=1 [V]\) on the inner boundary, and \(\Phi = 0 [V]\) on the outer boundary.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) | \(\|e\|_{H^1}\) | \(\alpha_{H^1}\) |
|---|---|---|---|---|---|---|---|
| 1 | 15 | 784 | 840 | 2.14e-03 | - | 1.08e-01 | - |
| 1 | 16 | 900 | 960 | 1.86e-03 | 2.00 | 1.01e-01 | 1.00 |
| 1 | 17 | 1024 | 1088 | 1.64e-03 | 2.00 | 9.43e-02 | 1.00 |
| 1 | 18 | 1156 | 1224 | 1.45e-03 | 2.00 | 8.87e-02 | 1.00 |
| 2 | 15 | 784 | 3248 | 3.80e-06 | - | 6.70e-04 | - |
| 2 | 16 | 900 | 3720 | 3.08e-06 | 3.07 | 5.83e-04 | 2.01 |
| 2 | 17 | 1024 | 4224 | 2.53e-06 | 3.06 | 5.12e-04 | 2.01 |
| 2 | 18 | 1156 | 4760 | 2.10e-06 | 3.05 | 4.53e-04 | 2.01 |
| 3 | 15 | 784 | 7224 | 7.00e-07 | - | 8.11e-05 | - |
| 3 | 16 | 900 | 8280 | 5.31e-07 | 4.00 | 6.59e-05 | 3.00 |
| 3 | 17 | 1024 | 9408 | 4.10e-07 | 4.00 | 5.43e-05 | 3.00 |
| 3 | 18 | 1156 | 10608 | 3.22e-07 | 4.00 | 4.53e-05 | 3.00 |
The following notations were used in the table header:
There is another interesting observation worth mentioning. The first convergence table above suggests that the \(L^2\) error norms are approximately the same for different degrees of interpolating Lagrange polynomials, \(p\), if the mesh refinement parameter, \(r\), is fixed. The normal behavior of the \( L^2 \) error norm, i.e., the higher the polynomial degree, the lower the error norm, is recovered by application of the exact boundary conditions, see the second convergence table above.
The figure below illustrates the electric potential calculated for \(a=0.5 [m] \), \(b=1.0 [m] \), and \(\Phi_0 = 1 [V] \).
The results of the three-dimensional version of the experiment are similar to those of the two-dimensional version of the experiment described above. The straightforward application of the Dirichlet boundary condition, i.e., running the executable file cbnd-shell, yields the following convergence table.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) | \(\|e\|_{H^1}\) | \(\alpha_{H^1}\) |
|---|---|---|---|---|---|---|---|
| 1 | 9 | 3072 | 3474 | 1.22e-02 | - | 3.37e-01 | - |
| 1 | 10 | 4374 | 4880 | 9.63e-03 | 1.98 | 3.00e-01 | 1.00 |
| 1 | 11 | 6000 | 6622 | 7.81e-03 | 1.99 | 2.70e-01 | 1.00 |
| 1 | 12 | 7986 | 8736 | 6.46e-03 | 1.99 | 2.45e-01 | 1.00 |
| 2 | 9 | 3072 | 26146 | 1.46e-02 | - | 8.21e-02 | - |
| 2 | 10 | 4374 | 36974 | 1.16e-02 | 1.98 | 6.85e-02 | 1.55 |
| 2 | 11 | 6000 | 50442 | 9.37e-03 | 1.99 | 5.82e-02 | 1.54 |
| 2 | 12 | 7986 | 66838 | 7.76e-03 | 1.99 | 5.03e-02 | 1.54 |
| 3 | 9 | 3072 | 86450 | 1.46e-02 | - | 7.90e-02 | - |
| 3 | 10 | 4374 | 122528 | 1.16e-02 | 1.98 | 6.58e-02 | 1.55 |
| 3 | 11 | 6000 | 167462 | 9.38e-03 | 1.99 | 5.60e-02 | 1.54 |
| 3 | 12 | 7986 | 222224 | 7.76e-03 | 1.99 | 4.83e-02 | 1.54 |
The following notations were used in the table header:
Here again, we can observe a reduction of the convergence orders in the comparison to the convergence orders observed in the Method of manufactured solutions (mms/) numerical experiment. The convergence orders can be recovered by application of the exact boundary conditions. That is, by running the executable file cbnd-shell-exact. The results of such a run are listed below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) | \(\|e\|_{H^1}\) | \(\alpha_{H^1}\) |
|---|---|---|---|---|---|---|---|
| 1 | 9 | 3072 | 3474 | 1.22e-02 | - | 3.37e-01 | - |
| 1 | 10 | 4374 | 4880 | 9.63e-03 | 1.98 | 3.00e-01 | 1.00 |
| 1 | 11 | 6000 | 6622 | 7.81e-03 | 1.99 | 2.70e-01 | 1.00 |
| 1 | 12 | 7986 | 8736 | 6.46e-03 | 1.99 | 2.45e-01 | 1.00 |
| 2 | 9 | 3072 | 26146 | 7.92e-05 | - | 7.76e-03 | - |
| 2 | 10 | 4374 | 36974 | 5.50e-05 | 3.10 | 6.12e-03 | 2.01 |
| 2 | 11 | 6000 | 50442 | 3.97e-05 | 3.09 | 4.95e-03 | 2.01 |
| 2 | 12 | 7986 | 66838 | 2.96e-05 | 3.07 | 4.09e-03 | 2.01 |
| 3 | 9 | 3072 | 86450 | 1.47e-05 | - | 1.04e-03 | - |
| 3 | 10 | 4374 | 122528 | 9.18e-06 | 3.99 | 7.29e-04 | 3.00 |
| 3 | 11 | 6000 | 167462 | 6.03e-06 | 3.99 | 5.32e-04 | 3.00 |
| 3 | 12 | 7986 | 222224 | 4.12e-06 | 3.99 | 3.99e-04 | 3.00 |
The following notations were used in the table header:
We can conclude the following:
If a boundary of a problem domain is polygonal, it can be approximated by the mesh exactly. In such cases, assuming that the solution is expected to be a smooth function without singularities, no question arises on how to allocate the computer resources: allocate them for increasing the degree of the interpolating polynomials, rather then on the increasing the amount of cells in the mesh. This conclusion is a result of contemplating the convergence tables of the Method of manufactured solution (mms/) numerical experiment.
If the boundaries of the problem domain are curved, the choice of how to allocate the computer resources is far from being obvious. In the two examples above increase of the degree of interpolating Lagrange polynomials has not yielded a decrease of the \(L^2\) error norm, see the first and the third convergence tables above. Therefore, increasing the amount of cells seems to be the optimal choice when allocating computer resources. We can conclude that a proposition to increase the degree of the interpolating polynomials must be treated with a healthy skepticism, unless the mesh approximates the boundaries of the problem domain exactly.
Note, in this particular experiment the degradation of the convergence rate can be avoided by invoking second-order mapping and attaching a spherical manifold to the mesh.
[1] Zienkiewicz, Olek C and Taylor, Robert Leroy and Zhu, Jian, "The finite element method: its basis and fundamentals", Elsevier, 2005.
1.9.1