Introduction

In the Method of manufactured solutions (mms/) numerical experiment, the volume free-charge density, $\rho_f$, on the right-hand side of the partial differential equation is a continuous function of spatial coordinates that fills all the space in the problem domain. In this numerical experiment we will test if StaticScalarSolver::Solver can handle volume free-charge densities with discontinuity.

We will do the test by the method of exact solutions. That is, we will solve two textbook problems in which $\rho_f$ is confined to a finite region of a problem domain and compare the numerical results with the closed-form analytical solutions.

Implementation

Two-dimensional problem

The following textbook problem is solved in the two-dimensional version of the experiment. (Section 6.1.11. Volume charge with a cylindrical interface)

Free charge is uniformly distributed within an infinity long cylinder of a radius $a$. This cylinder is enclosed in an infinitely long conducting tube with an inner radius of $b$. The cylinder and the tube are coaxial. The tube is grounded. All materials inside the conducting tube have permittivity of $\epsilon_0$. The figure below illustrates this configuration.

We would like to calculate the electric potential.

The solution to this problem reads

\begin{equation} \Phi = \left\{ \begin{aligned} & \frac{\rho_f a^2}{2 \epsilon_0} \ln{\frac{b}{r}} &\text{if } \text{ } & a \le r \le b \\ & \frac{\rho_f a^2}{4 \epsilon_0} \Bigg[1 + 2\ln{\frac{b}{a}} - \frac{r^2}{a^2} \Bigg] &\text{if } \text{ } & r \le a, \\ \end{aligned} \right. \end{equation}

\begin{equation} \vec{\nabla}\Phi = \left\{ \begin{aligned} & - \frac{\rho_f a^2}{2 \epsilon_0}\frac{1}{r} \hat{r} &\text{if } \text{ } & a \le r \le b \\ & - \frac{\rho_f}{2\epsilon_0} r \hat{r} &\text{if } \text{ } & r \le a. \\ \end{aligned} \right. \end{equation}

As soon as this problem exhibits translation symmetry, we will solve it numerically in a two-dimensional planar domain. We remove the conducting tube from the problem domain and model it by setting up a Dirichlet boundary condition $\Phi = 0$ on the outer boundary of the problem domain.

Three-dimensional problem

The following textbook problem is solved in the three-dimensional version of the experiment. (Section 6.1.12. Volume charge with a spherical interface)

Free charge is uniformly distributed within a spherical ball of a radius $a$. The ball is enclosed in a conducting spherical shell with an inner radius of $b$. The ball and the shell are concentric. The shell is grounded. All materials inside the conducting spherical shell have permittivity of $\epsilon_0$. The figure below illustrates this configuration.

We would like to calculate the electric potential.

The solution to this problem reads

\begin{equation} \Phi = \left\{ \begin{aligned} &\frac{\rho_f a^3}{3 \epsilon_0}\Bigg[\frac{1}{r} - \frac{1}{b}\Bigg] &\text{if } \text{ } & a \le r \le b \\ &\frac{\rho_f}{6 \epsilon_0} \Bigg[a^2 + 2a^3 \Bigg(\frac{1}{a} - \frac{1}{b}\Bigg) - r^2 \Bigg] &\text{if } \text{ } & r \le a, \\ \end{aligned} \right. \end{equation}

\begin{equation} \vec{\nabla}\Phi = \left\{ \begin{aligned} &- \frac{\rho_f a^3}{3 \epsilon_0} \frac{1}{r^2}\hat{r} &\text{if } \text{ } & a \le r \le b \\ &- \frac{\rho_f}{3\epsilon_0} r \hat{r} &\text{if } \text{ } & r \le a. \\ \end{aligned} \right. \end{equation}

We will solve this problem numerically in a three-dimensional domain. We remove the conducting shell from the problem domain and model it by setting up a Dirichlet boundary condition $\Phi = 0$ on the outer boundary of the problem domain.

The boundary value problem

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)=\rho_f & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &\Phi = 0 & \text{on} & \Gamma_{D1} & \text{(ii)},\\ \end{array} \end{equation}

where

\begin{equation} \rho_f = \left\{ \begin{aligned} & \rho_0 &\text{if } \text{ } & r < a \\ & 0 &\text{if } \text{ } & r > a. \end{aligned} \right. \end{equation}

This boundary value problem is a special case of the general static scalar boundary value problem.

The program

This experiment is implemented in accordance with the base code structure. This time we use second-order mapping and attach a spherical manifold to all faces of the mesh. That is, we declare the spherical manifold in the SolverRHO template,

template<int dim>
class SolverRHO
  : public SettingsRHO
  , public Solver<dim>
{
...
  SphericalManifold<dim> sphere;
...
};

and attach it to the mesh as the following:

template<int dim>
void
SolverRHO<dim>::mark_materials()
{
  Solver<dim>::triangulation.reset_all_manifolds();
 
  for (auto cell : Solver<dim>::triangulation.active_cell_iterators()) {
    
        ...
 
    for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; f++) {
      double dif_norm = 0.0;
      for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; v++) 
        dif_norm += std::abs(cell->face(f)->vertex(0).norm() -
                               cell->face(f)->vertex(v).norm());
 
      if ((dif_norm < eps) && (cell->center().norm() > d1))
        cell->face(f)->set_all_manifold_ids(1);
    }
  }
 
  Solver<dim>::triangulation.set_manifold(1, sphere);
}

That is, the spherical manifold is attached to a face if all the vertices of the face are equidistant from the origin and the center of the cell the face belongs to is outside a circle (ball) of a radius of $d1$ (the half-side of the square (cube) in the middle of the mesh). Then we pass 2 as the second parameter to the constructor of SolverRHO class template to enable second-order mapping,

class BatchRHO : public SettingsRHO
{
public:
  void run()
  {
...
          SolverRHO<DIMENSION__> problem(
            p,
            2,
            r,
            dir + "solution_p" + std::to_string(p) + "_r" + std::to_string(r));
...
    };
};

The build process generates two executable files: rho-circle and rho-sphere. That is, one executable file for the two-dimensional experiment, rho-circle, and one executable file for the three-dimensional experiment, rho-sphere. To rebuild them change into rho/build/Release directory and execute the following:

user@computer .../rho/build/Release$ ./clean

user@computer .../rho/build/Release$ ./build

Then all executable files must be executed again. This can be done by changing into the rho/bin/Release directory and executing run-all script there,

user@computer .../rho/build/Release$ cd ../../bin/Release

user@computer .../rho/bin/Release$ ./run-all

This will generate various files in the rho/bin/Release/Data directory. Among the generated files there are vtu files that can be viewed with a help of ParaView software package of the Kitware, Inc. The data files in tex and txt format contain the convergence tables. The rho/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 rho/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into rho/gmsh directory and executing the following:

user@computer .../rho/gmsh$ ./clean

user@computer .../rho/gmsh$ ./build

This will generate a set of globally refined meshes in rho/gmsh/data.

The SettingsRHO 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. The last can be useful when debugging.

Simulation results

All calculation have been done with the following settings, $\rho_0=\epsilon_0$, $a=0.5[m]$ and $b=1.0[m]$, see SettingsRHO.

Electric potential

Two-dimensional experiment

The plot below illustrates the numerically calculated electric potential.

The table below illustrates the error norms and the convergence orders.

p	r	cells	dofs	$\\|e\\|_{L^2}$	$\alpha_{L^2}$	$\\|e\\|_{H^1}$	$\alpha_{H^1}$
1	5	320	337	5.72e-04	-	1.80e-02	-
1	6	500	521	3.70e-04	1.96	1.45e-02	0.96
1	7	720	745	2.58e-04	1.97	1.22e-02	0.97
1	8	980	1009	1.91e-04	1.97	1.05e-02	0.98
2	5	320	1313	1.33e-05	-	6.93e-04	-
2	6	500	2041	6.88e-06	2.96	4.47e-04	1.97
2	7	720	2929	4.00e-06	2.98	3.11e-04	1.98
2	8	980	3977	2.52e-06	2.98	2.29e-04	1.98
3	5	320	2929	8.30e-07	-	4.07e-05	-
3	6	500	4561	3.42e-07	3.97	2.11e-05	2.95
3	7	720	6553	1.65e-07	3.98	1.23e-05	2.97
3	8	980	8905	8.95e-08	3.99	7.74e-06	2.98

The following notations were used in the table header:

p - the degree of the interpolating Lagrange polynomials that constitute the shape functions.
r - the number of nodes on transfinite lines, see discussion in here.
cells - the total amount of active cells.
dofs - the amount of degrees of freedom in the active cells.
$\|e\|_{L^2}$ - the $L^2$ error norm.
$\|e\|_{H^1}$ - the $H^1$ error norm.
$\alpha_{L^2}$ - the order of convergence of the $L^2$ error norm.
$\alpha_{H^1}$ - the order of convergence of the $H^1$ error norm.

Three-dimensional experiment

The plot below illustrates a numerically calculated electric potential.

The table below illustrates the error norms and the convergence orders generated by the three-dimensional version of the numerical experiment.

p	r	cells	dofs	$\\|e\\|_{L^2}$	$\alpha_{L^2}$	$\\|e\\|_{H^1}$	$\alpha_{H^1}$
1	5	832	909	7.63e-04	-	2.08e-02	-
1	6	1625	1736	5.02e-04	1.88	1.68e-02	0.94
1	7	2808	2959	3.59e-04	1.84	1.43e-02	0.91
1	8	4459	4656	2.66e-04	1.94	1.23e-02	0.98
2	5	832	6905	2.52e-05	-	1.25e-03	-
2	6	1625	13371	1.28e-05	3.03	8.09e-04	1.96
2	7	2808	22981	7.40e-06	3.02	5.64e-04	1.97
2	8	4459	36359	4.65e-06	3.01	4.16e-04	1.98
3	5	832	22981	9.01e-06	-	1.44e-04	-
3	6	1625	44656	3.73e-06	3.96	7.45e-05	2.96
3	7	2808	76915	1.80e-06	3.98	4.30e-05	3.01
3	8	4459	121864	9.78e-07	3.98	2.73e-05	2.96