 |
Sapphire++
v1.3.0-121-g2682132 |
|
Sapphire++
v1.3.0-121-g2682132 |
Working with discontinuous Galerkin (DG) finite element methods can be counterintuitive when it comes to visualizing and interpreting the results. As the name suggests, the solution is multivalued and discontinuous at cell edges, a feature which contrasts the nature of Finite Difference (FD) or Finite Volume (FV) methods. These methods by definition only have one value per node/cell. In this section, we want to provide you some tips, tricks and lessons we had to learn on the way, that might help with visualizing and interpreting the output of Sapphire++.
In finite element (FE) methods, differential equation are solved by finding weak solutions to the problem. This means, that the numeric solution \(f_h(x)\) does not need to converge to the full solution \(f\) point wise, but only in the integral sense. We express this in the following way:
\[ (f, v) = (f_h, v) \quad \forall v \in V_h \,, \]
where \(v\) is a test function and \((f, v) = \int f(x) v(x) dx\) is the inner product of the function space \(V_h\). Using that \(V_h\) is spanned by a set of basis functions \(\phi_i(x)\), Galerkin methods use them to project the solution \(f_h\) into the function space \(V_h\):
\[ f_h(x) = \sum_i f_i \phi_i(x) \,. \]
The numeric solution is now represented by the coefficients \(f_i\). We say that \(f_h(x)\) is the projection of \(f(x)\) onto the FE space. In the picture below we show the difference between a projection \(f_h\) and an interpolation. In the interpolation, the function has the exact value at the nodes of the mesh, using a linear interpolation between the nodes. This is not the case for the projection. In fact, using a discontinuous Galerkin method, the basis functions are only defined on one cell each. Therefore, the projection is discontinuous at the cell boundaries.
Running Sapphire++, we observed that at later time steps \(t^n\), these discontinuities might be more pronounced in the numerical solution \(f_h^n\), than those of the projection at the given time step \(f_h(t^n)\). This becomes especially noticeable, if the analytic solution would be zero at some time, \(f(t^n, x) = 0\). An example is the \(f_{110}\) mode in the convergence-study example.
These jumps are a feature of the DG method and are not a sign of instability. In fact, error estimates for DG methods include estimated for the jump terms, see [5]. A typical example is the \(L^2\) error estimate of a DG method is,
\[ \|u^N - u_h^N\|_{L^2} + \left( \sum_m^{N-1} \delta t \sum_F \int_F \frac{1}{2} |\beta \cdot n_F| [[u^m - u_h^m]]^2 \right)^{1/2} \lesssim C (\Delta x)^a + D (\Delta t)^b \,. \]
The second term on the left-hand side provides an estimate of the jumps at the cell boundaries, \([[u]] = \sum_j (u(x^+_{j+\frac{1}{2}}) - u(x^-_{j+\frac{1}{2}}))\). \(a\), \(b\), \(C\), and \(D\) are constants that bound the error. Concluding, we can see that the discontinuities are a feature of the DG method, that we can expect to be on the order of numerical errors.
In Sapphire++, we use polynomial basis functions \(\phi_i\) for the FE representation. More precisely, we Lagrange polynomials of degree \(k\) on each cell. Representing these in a visualization programme, like ParaView and VisIt, turn out to be challenging. Most programmes expect a point wise representation of the solution, while the FE representation is given by the coefficients \(f_i\). Hence, when saving the results deal.II "converts" it to point wise data, by giving the values on each node of the FE mesh. In the special case of first order polynomials \(k=1\), this results in an exact representation of the numerical solution, as most visualisation programmes will linearly interpolate between the nodes. For higher order polynomials, deal.II provides the following workaround: it adds additional nodes to the mesh, representing the polynomials inside the cells. There are some important points to note:
This figure shows again the \(f_{110}\) mode of the convergence-study example, but with a higher polynomial degree of \(k=2\). One can see that the discontinuities much less pronounced than in the \(k=1\) case (figure in previous section).
In the VFP module in Sapphire++, we have to consider some additional points when visualizing and interpreting the results:
Here we want to share some tips and tricks we learned when using ParaView to visualize the results of Sapphire++:
To extract a 1D line from 2D/3D data, use the PlotOverLine feature. The best way to use this feature is to select the Sample At Cell Boundaries option to get the full discontinuous representation of the data. If you prefer a "smooth" result, the Sample At Segment Centers option is a good choice. (For \(k=1\), this equals the cell average.)
Unfortunately, the PlotOverLine feature has a few quirks we want to mention here:
1.16.1