Dear community,
I have a general question for the adaptive finite element method. I’m new to the finite element method and numerical solvers so I’m really sorry if this is stupid. The short version is:
When testing an error estimator on a problem without an exact solution, some literature uses the numerical solution from the final adaptive refinement step as a reference solution (u_{ref}). To compute the error \|u_{ref}-u_h\| of another numerical solution (u_h), they interpolate u_h into u_{ref}'s function space. Is this procedure theoretically sound? Doesn’t it assume the adaptively refined mesh is already reliable, creating potential circular numerical reasoning?
The detailed context with a concrete demo is as follows:
Suppose we want to solve a problem on a unit cube starting with initial mesh:
We perform 7 rounds of uniform refinement to obtain a series of uniformly refined meshes and solutions. The last uniformly-refined mesh:
Say now we use AFEM (starting from the same initial mesh) following the classic recipe SOLVE-ESTIMATE-MARK-REFINE (in this demo, we are actually keeping mark the cylindrical region \{(x, y, z): (x-1)^2+(y-1)^2<0.1^2\}) and we get the last adaptively-refined mesh:
We set the reference solution u_{ref}:=u^{(11)}_{adapt} as the last adaptively-solved solution. To compute H(curl) error of the numerical solution u_h which is from uniform mesh hierarchy or adapt mesh hierarchy, I hope that we can do the following:
from dolfinx.fem import (
Function,
form, assemble_scalar, create_interpolation_data,
)
from ufl import (
inner, curl, dx,
)
from mpi4py import MPI
import numpy as np
def calc_hcurl_error(u, uh, ref=False):
if ref:
domain_fine = u.function_space.mesh
uh.x.scatter_forward()
padding = 1e-16
fine_mesh_cell_map = domain_fine.topology.index_map(domain_fine.topology.dim)
num_cells_on_proc = fine_mesh_cell_map.size_local + fine_mesh_cell_map.num_ghosts
cells = np.arange(num_cells_on_proc, dtype=np.int32)
interpolation_data = create_interpolation_data(u.function_space, uh.function_space, cells, padding=padding)
uh_interp = Function(u.function_space)
uh_interp.interpolate_nonmatching(uh, cells, interpolation_data)
uh = uh_interp
error_form = form(inner(curl(u - uh), curl(u - uh))*dx + inner(u - uh, u - uh)*dx)
error = np.sqrt(uh.function_space.mesh.comm.allreduce(assemble_scalar(error_form), MPI.SUM))
return error
My first question is, is it valid to use the final adaptive solution as the reference? For instance, comparing errors ||u_{adapt}^{(10)} - u_{ref}|| and ||u_{uniform}^{(7)} - u_{ref}||, where u_{ref} is u_{adapt}^{(11)}, seems problematic. Since the meshes for u_{adapt}^{(10)} and u_{adapt}^{(11)} are very similar, their difference might be artificially small, potentially biasing the comparison in favor of the adaptive method. Doesn’t this resemble a cyclic argument?
In my own opinion, a more mathematically rigorous, but computationally much heavier alternative could be ‘merging’ two meshes, or the common refinement, i.e. find a mesh \mathcal{M}_f which is strictly finer than both the reference mesh \mathcal{M}_{ref} and the mesh of numerical solution \mathcal{M}_{h}. So we could do prolongations (in dolfinx I think it could be done implicitly by interpolate_nonmatching) from \mathcal{M}_{ref} to \mathcal{M}_f and from \mathcal{M}_{h} to \mathcal{M}_{f}. However, in practice, I believe that exactly merging two arbitrary, non-nested meshes (with zero tolerance) is computationally prohibitive.
So besides interpolation to a reference adaptive mesh, are there other procedures considered best practice for the referenced error computation when an exact solution is unavailable? Are there any built-in functions in Dolfinx or Netgen that support computing the intersection or common refinement of two existing volume meshes? I’m currently trying Salome to merge intersecting meshes Building Compound Meshes — Mesh 9.15.0 documentation.
Any comments, ideas, or guidance from your experience would be greatly appreciated. Thank you for your time and patience in reading through this！