In section 2.2 of reference [1], the authors say that a way to calculate the rate of convergence of a solution is to calculate a sequence of ratios as follows:
where u is the numerical solution, and h is mesh size. These ratios should be about the same in the asymptotic region. For a smooth error, this rate should be the same rate of convergence as when compared to the actual solution (if known).
How do I go about properly calculating this type of ratio using FEniCS, given how there are subtleties when it comes to choosing which points the numerical solution is evaluated at? My first thought is to first interpolate the solution and use that at a few select points within the mesh, but I just wanted to make sure there is some standard way out there to tackle this first.
I understand that I can use the known exact solution in errornorm() to calculate various types of error but I am curious about this method. Any advice is appreciated! I am relatively new to FEniCS.
you have already noticed that computing the ratios can be difficult. Actually note that there should be some kind of norm taken, i.e., you want to compare the ratios of
| u_h - u_h/2 | / | u_h/2 - u_h/4 |
and so on, where | * | denotes the norm in whose convergence you are interested in.
In order to calculate this properly (or “exactly”) you need a mesh hierarchy. I explain this using the unit square and quadriliterals as an example, and I assume that the discretization is uniform, the generalization of the idea to complicated meshes and other elements is in theory straightforward, in practice it may not.
So assume you have the mesh hierarchy from a uniform refinement, i.e., you split every quadriliteral into 4 others evenly. Then, every grid point of your coarse mesh (corresponding to size h) is also a grid point of your fine mesh (corresponding to size h/2). Now, if you interpolate the solution on the coarse mesh onto the fine mesh, you get exactly the same solution, atleast for Lagrangian elements, as the interpolation only depends on function evaluations (at the mesh points and the dofs - for higher order elements). Therefore, you can now take the difference of these solutions on the fine mesh and everything is well-defined.
However, as soon as you don’t have this uniform refinement type mesh hierarchy, the things are not so nice and in theory all kinds of unwanted effects can happen (you can think of easy examples in 1D for piecewise linear functions). But if your mesh is sufficiently fine it should not matter that much.
So if possible, I would try to generate the meshes from uniform refinement, but you can also just proceed with finer meshes that are not necessarily obtained by splitting the elements. Still, I would always interpolate (and not project) the solution on the coarse mesh onto the fine mesh and compare the solutions there as you have “more information” available there and in general this should work nicely.
Thank you for your in-depth response! I am trying to apply your response to a case involving Nedelec edge elements of the first kind on tetrahedra in a unit cube mesh.
However, I am having trouble visualizing uniform refinement in this scenario. Will comparing interpolations between fine and coarse meshes still make sense?
For instance, my coarse mesh could be
mesh = UnitCubeMesh(h, h, h)
and my fine mesh
mesh = UnitCubeMesh(2*h, 2*h, 2*h)
To my understanding, h and 2*h above relates to a mesh resolution of sorts. Will the above definitions of ‘mesh’ necessarily translate to uniform refinement? Apologies if this is a trivial question.
You’re welcome.
I actually have no experience with Nedelec elements, but maybe someone else can help you with that.
For your question regarding the mesh: As far as I know, the unit* meshes generate structured grids. Your idea works just fine. However, I would not specify the mesh using “h”, which is in my understanding a length scale, as the input parameters for unit* meshes should actually be integers that describe how many elements there will be along a certain axis. Let me describe this in 2D for the UnitSquareMesh, it extends naturally to 3D. The command
UnitSquareMesh(n_x, n_y)
generates a mesh of the UnitSquare by meshing it with quadrilaterals. It uses n_x quads in x direction and n_y quads in y direction. Then, the quads are split in half such that we get two triangles from each quad. Therefore, using
UnitSquareMesh(2*n_x, 2*n_y)
gives us a uniform refinement of the mesh. Now, the same applies in 3D with the UnitCubeMesh, where every hexahedron is split into 6 tetrahedrons.