 |
Logbook
(07-04-2025)
Static problems
|
|
|
Logbook
(07-04-2025)
Static problems
|
|
Classes:
1) BatchMMS
2) SolverMMS
3) ExactSolutionMMS_PHI
4) ExactSolutionMMS_E
5) ExactSolutionMMS_D
6) SettingsMMS
Files:
1) mms/src/main.cpp
2) mms/include/solver.hpp
3) mms/src/solver.cpp
4) mms/include/exact_solution.hpp
5) mms/src/exact_solution.cpp
6) mms/src/static_scalar_input.cpp
7) mms/include/settings.hpp
The method of manufactured solutions (MMS) is a method of verification of scientific software. It verifies if a code indeed solves a given boundary value problem. It does not, however, validate the boundary value problem. Neither does it validate the formulation of the problem. Simply put, the MMS makes sure that the code does not contain bugs. The purpose of this numerical experiment is to verify the code of the following class templates:
The first class template implements a solver that solves the static scalar boundary value problem which is based in the div-grad partial differential equation. The other two class templates convert the calculated electric potential, \(\Phi\), i.e., the output of StaticScalarSolver::Solver, into electric field, \(\vec{E}\), and displacement, \(\vec{D}\).
In the case of the StaticScalarSolver::Solver class template, the MMS 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 \(\Phi_M\). Second, we differentiate the manufactured solution several times to get closed-form expressions of the right-hand sides of the boundary value problem, i.e., \(\rho_f\), \(\eta\), and \(\sigma\). We choose closed-form expressions for \(\epsilon\) and \(\gamma\) arbitrary. Third, we obtain a numerical solution, \(\Phi\), by feeding to the solver the closed-form expressions for \(\rho_f\), \(\eta\), \(\sigma\), \(\epsilon\), and \(\gamma\). Fourth, we calculate an error as a difference between the closed-form manufactured solution and the numerical solution:
\[ e = \Phi_M-\Phi. \]
Finally, we evaluate the \(L^2\) and \(H^1\) norms of the error. The \(L^2\) norm is evaluated as the following
\[ ||e||_{L^2} = \sqrt{\iiint_{\Omega} e^2 dV}. \]
The \(H^1\) seminorm is, essentially, the \(L^2\) norm of the gradient of the error:
\[ ||e||_{H^1} = \sqrt{\iiint_{\Omega} \big( \vec{\nabla}e \big)^2 dV}. \]
The last two equations are valid for three-dimensional problems. The triple volume integrals must be replaced with double integrals over the surface of the problem domain in two-dimensional problems. We repeat the calculations of the error norms on a set of meshes that are refined versions of a base coarse mesh. The \(L^2\) and \(H^1\) error norms are supposed to converge to zero at rates predicted by the theory as the size of the mesh cells decreases. If they do converge to zero at predicted rates, we conclude that the amount of bugs in the code is acceptable.
In theory, see the Lecture 3.91 and Lecture 3.95 by Wolfgang Bangerth, both error norms should converge to zero at the following convergence rate:
\[ ||e|| = c h^{\alpha}, \]
where \(h\) is the size of the mesh cells, \(c\) is a constant, and \(\alpha\) is the convergence order. The theory predicts the following convergence orders:
\[ \alpha_{L^2} \le p + 1 \]
for the \(L^2\) error norm, and
\[ \alpha_{H^1} \le p \]
for the \(H^1\) error norm. The parameter \(p\) in the last two equations is a degree of the Lagrange interpolating polynomials that constitute the FE_Q finite elements.
We will verify the StaticScalarSolver::ProjectPHItoE and StaticScalarSolver::ProjectPHItoD class templates as the following. First, we will obtain manufactured solutions \(\vec{E}_{M}\) and \(\vec{D}_M\), by differentiating \(\Phi_M\), i.e.,
\[ \vec{E}_M = - \vec{\nabla} \Phi_M \]
and
\[ \vec{D}_M = - \epsilon \vec{\nabla} \Phi_M. \]
Second, we will feed to StaticScalarSolver::ProjectPHItoE and StaticScalarSolver::ProjectPHItoD the numerical solution \(\Phi\) and the manufactured solutions \(\vec{E}_M\) and \(\vec{D}_M\) to obtain two \(L^2\) error norms, one for the electric field and one for the displacement. If the error norms will converge to zero at any rate as we refine the mesh, we will conclude that the amount of bugs in the code is acceptable.
The MMS has a drawback: it cannot be used for verification of all aspects of a boundary value problem within one numerical experiment. For example, it is quite difficult to manufacture a solution that would incorporate a non-trivial configurations of the volume free charge, \(\rho_f\), the surface free charge, \(\kappa_f\), and discontinuity of the permittivity, \(\epsilon\), on interfaces between dielectrics. One way to approach this problem is to split verification of all aspects of the boundary value problem between a number of numerical experiments. To do so, we set \(\kappa_f\) to zero in this numerical experiment and assume that the permittivity, \(\epsilon\), is a smooth function of the spatial coordinates. The surface free charge, \(\kappa_f\), is tested in the Surface charge (sch/) numerical experiment. The discontinuity of the permittivity, \(\epsilon\), is put to a test in the Interface between dielectrics (int/) numerical experiment.
We chose (somewhat arbitrary) the following manufactured solution
\[ \Phi_M(x, y) = \cos{(kx)} + \cos{(ky)}, \]
and the coefficients on the left-hand sides of the boundary value problem,
\[ \epsilon = \epsilon_0 (x^2 y^2 + 1), \]
\[ \gamma = \epsilon \big( \sqrt{x^2 + y^2} + 2 \big). \]
The parameter \(k\) can be adjusted in the code, i.e., in class SettingsMMS. The gradient of the manufactured solution and the right-hand side of the partial differential equation, \(\rho_f\), is obtained by straightforward differentiation of the equations above:
\[ \vec{\nabla} \Phi_M(x, y) = -k \sin{(kx)} \hat{i} -k \sin{(ky)}\hat{j}, \]
\[ \rho_f = -\vec{\nabla} \cdot \bigg(\epsilon \vec{\nabla} \Phi_M \bigg) =\epsilon_0 k \bigg\{ 2xy^2 \sin{(kx)} + 2x^2y\sin{(ky)} + k(x^2y^2+1)\big(\cos{(kx)} + \cos{(ky)} \big)\bigg\}. \]
The problem domain is a square with a side of two meters or a circle with a diameter of two meters. Both, the square and the circle are centered at the origin as shown in the figure below. Both problem domains are assumed to be planar, i.e., they are the result of exploiting a translation symmetry in the initial three-dimensional problem domain. 
The Dirichlet and Robin boundary conditions are applied on non-intersecting segments of the boundary. The boundary IDs are assigned in accordance with the boundary ID convention in the respective geo files. The square.geo file, for example, defines the boundary IDs as the following:
The simulation program simply loads these boundary IDs together with the mesh.
The right-hand side of the Dirichlet boundary condition, \(\eta\), is implemented simply by evaluating the manufactured solution on the boundary with boundary ID=1. That is,
\[ \eta = \Phi_M. \]
Similarly, the right-hand side of the Robin boundary condition is implemented by evaluating the following equation on the boundary with a boundary ID=2:
\[ \sigma = \epsilon \hat{n} \cdot \vec{\nabla} \Phi_M + \gamma \Phi_M. \]
A straightforward differentiation yields the manufactured solutions for the electric field and displacement:
\[ \vec{E}_M = - \vec{\nabla} \Phi_M(x, y) = k \sin{(kx)} \hat{i} + k \sin{(ky)} \hat{j}, \]
\[ \vec{D}_M = - \epsilon \vec{\nabla} \Phi_M(x, y) = \epsilon_0 (x^2 y^2 + 1) \big( k \sin{(kx)} \hat{i} + k \sin{(ky)} \hat{j} \big). \]
The solution in three dimensions is manufactured in a manner similar to that of the two-dimensional case. That is, we assume that the manufactured solution and the coefficients on the left-hand side of the boundary value problem have the following form
\[ \Phi_M(x, y, z) = \cos{(kx)} + \cos{(ky)} + \cos{(kz)}, \]
\[ \epsilon = \epsilon_0 (x^2 y^2 z^2 + 1), \]
\[ \gamma = \epsilon \big( \sqrt{x^2 + y^2 + z^2} + 2 \big). \]
Then the gradient of the manufactured solution and the right-hand side of the partial differential equation is calculated by differentiating the manufactured solution,
\[ \vec{\nabla} \Phi_M(x, y, z) = -k \sin{(kx)} \hat{i} -k \sin{(ky)}\hat{j} -k \sin{(kz)}\hat{k}, \]
\begin{equation} \begin{aligned} &\rho_f = -\vec{\nabla} \cdot \bigg(\epsilon \vec{\nabla} \Phi_M \bigg)=\\ &=\epsilon_0 k \bigg\{ 2xy^2z^2 \sin{(kx)} + 2x^2yz^2\sin{(ky)} + 2x^2y^2z\sin{(kz)}+k(x^2y^2z^2+1)\big(\cos{(kx)} + \cos{(ky)} + \cos{(kz)} \big)\bigg\}. \end{aligned} \end{equation}
The problem domain is a cube with a side of two meters or a sphere with a diameter of two meters. Both, the cube and the sphere are centered at the origin as shown in the figure below.
The boundaries of the domain are labeled in accordance with the boundary ID convention in the respective geo files. The cube.geo file, for example, defines the boundary IDs as the following:
The simulation program simply loads these boundary IDs together with the mesh.
The right-hand sides of the Dirichlet and the Robin boundary conditions are implemented by evaluating the following expressions on the corresponding boundaries
\[ \eta = \Phi_M, \]
\[ \sigma = \epsilon \hat{n} \cdot \vec{\nabla} \Phi_M + \gamma \Phi_M. \]
A straightforward differentiation yields the manufactured solutions for the electric field and displacement:
\[ \vec{E}_M = - \vec{\nabla} \Phi_M(x, y, z) = k \sin{(kx)} \hat{i} + k \sin{(ky)} \hat{j} + k \sin{(kz)} \hat{k}, \]
\[ \vec{D}_M = - \epsilon \vec{\nabla} \Phi_M(x, y, z) = \epsilon_0 (x^2 y^2 z^2 + 1) \big( k \sin{(kx)} \hat{i} + k \sin{(ky)} \hat{j} + k \sin{(kz)} \hat{k} \big). \]
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 \vec{\nabla} \Phi \big)= \rho_f & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &\Phi = \eta & \text{on} & \Gamma_{D1} & \text{(ii)},\\ \text{(n)}&\epsilon\hat{n}\cdot\vec{\nabla}\Phi+\gamma\Phi=\sigma&\text{on}&\Gamma_{R1}&\text{(iii)}.\\ \end{array} \end{equation}
This boundary value problem is a special case of the general static scalar boundary value problem. The boundary value problem is solved with a help of the StaticScalarSolver::Solver class template.
After the boundary value problem is solved, the numerical solution, \(\Phi\), is fed to the code implemented by the StaticScalarSolver::ProjectPHItoE class template which calculates the electric field as
\[ \vec{E} = - \vec{\nabla} \Phi. \]
Similarly, \(\Phi\), is fed to the code implemented by the StaticScalarSolver::ProjectPHItoD to compute the displacement as
\[ \vec{D} = - \epsilon \vec{\nabla} \Phi. \]
The MMS experiment is implemented in accordance with the base code structure. The build process generates four executable files: mms-square, mms-circle, mms-cube, and mms-sphere. To rebuild them change into mms/build/Release directory and execute the following:
Then all executable files must be executed again. This can be done by changing into the mms/bin/Release directory and executing run-all script there,
This will generate various files in the mms/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. Alternatively, the four executable files can be run one-by one.
Note that executable files require a set of meshes to be present in the mms/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into mms/gmsh directory and executing the following:
This will generate a set of refined meshes in mms/gmsh/data.
The SettingsMMS 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 vtk files next to the numerical solution. The last can be useful when debugging.
All the simulations described below have been done with a setting \(k=\pi\).
The four figures below illustrate the calculated electric potentials and list the corresponding convergence data.
mms-circle 
mms-cube 
mms-sphere 
The following notations are used in the headers of the tables that describe the convergence of the electric potential, \(\Phi\):
The values in the \(\alpha_{L^2}\) columns of the tables above approximately equal the values in the \(p\) columns plus one. This corresponds to the upper boundary of the theoretical prediction discussed in the introduction. Similarly, the values in the \(\alpha_{H^1}\) columns equal approximately the values in the \(p\) columns, what corresponds to the upper boundary of the theoretical predictions as well. We can conclude that the code of the solver implemented by the class template StaticScalarSolver::Solver contains an acceptable amount of bugs.
An analysis of the error norms listed in the tables above leads to a conclusion that it is always beneficial to allocate computer resources for increasing the degree of the interpolating polynomials, \(p\), rather than decreasing the size of the mesh cells, \(h\). This is, however, a wrong conclusion. In this particular experiment the convergence orders were at the upper boundary of the theoretical predictions, i.e.,
\[ \alpha_{L^2} \approx p + 1 \]
and
\[ \alpha_{H^1} \approx p \]
only because we have used the manufactured solution to calculate the boundary conditions at the exact locations of the support points. This is not possible in real-life situations. Consequently, the error norms in real-life simulations can converge slower than in this numerical experiment. Then the idea to decrease the size of the mesh cells, \(h\), can appear to be more attractive. The numerical experiment Effect of curved boundaries (cbnd/) allows to observe this effect. Furthermore, many problems in electromagnetics are formulated in unbounded domains. Direct implementation of an unbounded domain is impossible. For this reason, the boundary at infinity is replaced with a surface of finite dimensions. This replacement is also a kind of approximation that can decrease the convergence rate. Two methods of dealing with unbounded problem domains can be studied with a help of the numerical experiment Asymptotic boundary conditions (abc/).
The four figures below illustrate the calculated electric field, its magnitude, and the corresponding convergence data.
mms-square 
mms-circle 
mms-cube 
mms-sphere 
The following notations are used in the headers of the tables that describe the convergence of the electric field, \( \vec{E} \):
The four tables above suggest that the \(L^2\) error norm converges to zero as the size of mesh cells gets smaller. Therefore, we can conclude that the StaticScalarSolver::ProjectPHItoE class template contains an acceptable amount of bugs. It looks as the \(L^2\) convergence rate of the electric field, \(\vec{E}\), is similar to the \(H^1\) convergence rate of the electric potential, \(\Phi\): both of them approximately equal \(p\). The last is, most likely, due to the fact that the electric field and electric potential are related by the gradient.
The four figures below illustrate the calculated displacement, its magnitude, and the corresponding convergence data.
mms-square 
mms-circle 
mms-cube 
mms-sphere 
The following notations are used in the headers of the tables that describe the convergence of the displacement, \( \vec{D} \):
The four tables above suggest that the \(L^2\) error norm converges to zero as the size of mesh cells gets smaller. Therefore, we can conclude that the StaticScalarSolver::ProjectPHItoD class template contains an acceptable amount of bugs.
1.9.1 
