DWR Mesh-Refinement

Hi guys,

I’m looking through the “dwr_poisson.py” demo in Nate Sime’s dolfin_dg project dwr, and I’m a bit confused by one of the variables he uses. Particularly, lines 42 - 44 read

    j_h = Constant(10.0)*u*dx

lape = LinearAPosterioriEstimator(a, L, j_h, soln, bcs)


Presumably j_h is some sort of error measure - I’m not sure. Can someone clue me in?

1 Like

j_h is the functional of interest. I.e. it is the “goal” in the goal oriented dual weighted residual (DWR) adaptive refinement scheme. It is the quantity we want to optimise for error reduction. You can read more about that specific example here: https://fenicsproject.org/pub/course/lectures/2013-11-logg-fcc/lecture_11_error_control.pdf.

2 Likes

Thanks, Nate. One question though - from what I can tell, the functional of interest is determined by solving adjoint of the original BVP. However, what if we’re dealing with a problem that admits no analytical solution? Would it still make sense to use \int_{ \Omega} u \, dx as the functional of interest to minimize?

You don’t need to know any analytical solution to be able to use DWR. Sure, it might be beneficial to know the analytical value of J(u) to verify that the true error J(u) - J(u_h) and the DWR error estimator are similar. But even for this we frequently approximate J(u) by evaluating our quantity of interest with a FEM solution on a very fine mesh. From a mathematical point of view, concerning the regularity of the goal functional, there might be some restrictions. From a practical point of view, you can take almost any functional that you can come up with. E.g. we can also use point evaluations of the solution J(u) = u(x0) and DWR still tends to work, see e.g. https://arxiv.org/pdf/2404.01738.pdf Section 6.2.3 and Figure 5.

2 Likes

So I actually do have one more question about your implementation of DWR mesh-refinement. I was taking a closer look at the dwr_poisson.py script and I noticed something interesting -

The interior of the first for loop, where the mesh refinement is performed, reads as follows:

    V = FunctionSpace(mesh, 'CG', poly_o)
u, v = TrialFunction(V), TestFunction(V)
L = f*v*dx

soln = Function(V)
bcs = DirichletBC(V, Constant(0.0), 'on_boundary')
solve(a == L, soln, bcs)

dwr_errors.append(abs(assemble(soln*dx) - j))

dwr_dofs.append(V.dim())

j_h = Constant(10.0)*u*dx

lape = LinearAPosterioriEstimator(a, L, j_h, soln, bcs)
markers = lape.compute_cell_markers(FixedFractionMarker(frac=0.2))

mesh = refine(mesh, markers, redistribute=True)


My confusion lies in the fact that u is defined as a FEniCS TrialFunction, not the solution of the differential equation. That being said, when we evaluate the functional \int_{\Omega} 10 \,u \, dx, we’re not integrating the solution function, we’re just integrating an arbitrary element of the function space V. Am I correct or is there something I’m missing?

So I was doing a bit of experimenting with the dwr_poisson.py program and apparently DWR does not always work better than uniform mesh refinement. Based on the figure below, DWR converges faster than uniform refinement until ~120 DoFs have been used. After that, for a fixed number of DoFs, uniform refinement seems to perform much better than DWR. Is this simply a fundamental limitation of what we can expect from DWR?

The effectivities will depend on the scheme, regularity of the true solution, selection criteria, fraction of (de/)refinement, etc. One of the biggest aspects missing for exponential convergence is p-refinement. For problems with smooth solutions, adaptive h-refinement is not particularly great.

1 Like

On a semi-related note - is it possible to use point-evaluations in your DWR implementation? I tried using

j_h = u( (0.5, 0.5) )

as the functional of interest but that generated the error message,

Couldn't map 'v_1' to a float, returning ufl object without evaluation.


Are the only functionals supported by this framework UFL forms?

I’ve used a mollified delta function in the past which is somewhat useful.

Looks like I haven’t updated this example to dolfinx, but hopefully you get the idea. I can’t remember what tricks I used to get the linear estimator to work. So you may have to recast your problem as nonlinear. This shouldn’t be too difficult, and is equivalent to solving the linear problem.

In fact the linear estimator is just a recasting of a nonlinear estimator; dolfin_dg/dolfin_dg/dolfinx/dwr.py at ca7ea93216e7a871e30d7b07036919153f741029 · nate-sime/dolfin_dg · GitHub

I would love to spend more time on dolfin_dg, but it is not in the scope of any of my priorities.