Equation's strong residual is "big"

dolfinx
1

Hello, while fiddling with something loosely related, I noticed the following strange behaviour.
Consider the equation

-\Delta u = 1 - u^2 + f

where f and the Dirichlet boundary conditions are chosen so that the solution is, say, u_{ex} = \sin(12x) \cdot e^{-y^2} + \ln(1+z) [there is nothing special about this function].
If dolfin-x, I assemble the scalar residual = (div(grad(u)) + 1 - u**2 + f)**2*dx which is the equation’s strong residual in L2 norm.

Naively, I expect this quantity to be very small, but it turns out that it’s quite “big”: to give some sample output, on a unit cube with n subdivisions in each dimension, [after taking the square root]

n = 10; h = 0.173205;  ||L u_ex - f(u_ex)|| = 26.2191
n = 20; h = 0.0866025; ||L u_ex - f(u_ex)|| = 13.3575
n = 40; h = 0.0433013; ||L u_ex - f(u_ex)|| = 6.70934

Now, the first thing I notice is that the residual decreases linearly with h, but my question is on how to interpret that number: is 13 a lot? Is it small? Should I normalise it with something?

Changing the interpolation spaces (P2, P3, …) makes the residual smaller, but still not the “10^-many” I was expecting.
Changing the quadrature degree in UFL doesn’t seem to affect the computed residual at all.

Some context: I want to solve some nonlinear equations [like the one above] with a fixed-point method, using the equation’s strong residual as a stopping criterion.
In practice I see that, despite converging to the correct solution, the residual plateaus around the quantity computed with the exact solution [of course it does!], but then I am unable to set a meaningful tolerance for my stopping criterion.
To clarify the last paragraph: what I mean is that if I run the fixed-point procedure on the equation above with n = 10, the residuals I get for the various u_1, u_2, …, u_k make a sequence that converges to around 26. If I pick n = 20, I get another sequence that converges to around 13, and so on.

My underlying issue is that I want that residual to be very small, and I don’t understand why it isn’t :smiley:

Thanks for any help!

2

See for instance: Laplacian of solution extremely noisy

As you would expect from a variational solution, it does not satisfy the same continuity requirements as the strong formulation. I would for instance suggest you to read chapter 2 of the fenics tutorial

3

Using a straightforward “strong formulation residual” measured in the L_2 norm is typically not a good measure of the a posteriori error in the FE approximation of the true solution. Consider Section 4.1 of Graetsch & Bathe 2005.

Summarising: subject to the homogeneous boundary conditions of your example, it would be more informative to compute

\begin{align} \Vert u - u_h \Vert_E &= \sqrt{a(u - u_h, u - u_h)} \\ &\leq \left( \sum_{\kappa \in \mathcal{T}_h} c h^2_\kappa \Vert R \Vert^2_{L_2(\kappa)}\right)^\frac{1}{2} \end{align}

where R is the strong residual, h_\kappa is some geometric measure of an element in the mesh \kappa \in \mathcal{T}_h, a(\cdot, \cdot) is the appropriate bilinear component of your weak formulation and c is an unknown constant.

I wrote this pretty hastily, but cf. the following code (which could probably be made more efficient with some clever projection of the residual onto a DG0 space). We expect the error measured in the energy norm, and its estimate, to converge at an optimal rate of \mathcal{O}(h^p).

from dolfin import *
import numpy as np

n_eles = np.array([4, 8, 16])
a_prio_errs = np.zeros_like(n_eles, dtype=np.double)
a_post_errs = np.zeros_like(n_eles, dtype=np.double)
effectivities = np.zeros_like(n_eles, dtype=np.double)
R_l2s = np.zeros_like(n_eles, dtype=np.double)
h_sizes = 1.0/n_eles
p = 1

for i, n_ele in enumerate(n_eles):
	mesh = UnitSquareMesh(n_ele, n_ele)

	# FE problem
	V = FunctionSpace(mesh, "CG", p)
	u, v = TrialFunction(V), TestFunction(V)
	a = dot(grad(u), grad(v))*dx

	f = Expression("2*pi*pi*sin(pi*x[0])*sin(pi*x[1])", element=V.ufl_element())
	L = f*v*dx
	bc = DirichletBC(V, Constant(0.0), "on_boundary")

	u = Function(V)
	solve(a == L, u, bc)

	# `strong' residual
	R = div(grad(u)) + f

	# L2 norm of the residual
	R_l2s[i] = assemble(R**2*dx)**0.5

	# A posteriori unknown error energy norm
	# Unknown constant, but is approximately c = (1.0/(p+1))**2
	c = 1.0
	a_post_err = sum(c*(cell.circumradius())**2 * assemble_local(R**2*dx, cell) for cell in cells(mesh))**0.5
	a_post_errs[i] = a_post_err

	# A priori known error energy norm
	u_exact = Expression("sin(pi*x[0])*sin(pi*x[1])", degree=V.ufl_element().degree()+2, domain=mesh)
	a_prio_err = assemble(dot(grad(u_exact - u), grad(u_exact - u))*dx)**0.5
	a_prio_errs[i] = a_prio_err

	# Ratio of error estimate to true error
	effectivities[i] = a_post_err / a_prio_err


print("Residual L2 norms: ", R_l2s)
print("a posteriori error energy norms:", a_post_errs)
print("a priori error energy norms:", a_prio_errs)
print("effectivities:", effectivities)
print("a priori E norm convergence rates:", np.log(a_prio_errs[1:]/a_prio_errs[:-1])/np.log(h_sizes[1:]/h_sizes[:-1]))
print("a posteriori E norm convergence rates:", np.log(a_post_errs[1:]/a_post_errs[:-1])/np.log(h_sizes[1:]/h_sizes[:-1]))

with p=1 I have the following output:

Residual L2 norms:  [8.93204975 9.62080787 9.80649238]
a posteriori error energy norms: [1.57897824 0.85036731 0.43338983]
a priori error energy norms: [0.86352266 0.43581471 0.2180719 ]
effectivities: [1.82853133 1.95121294 1.9873713 ]
a priori E norm convergence rates: [0.98651913 0.99891105]
a posteriori E norm convergence rates: [0.89283324 0.97242084]

So we can see: The L_2 norm of the strong residual grows, telling us nothing about the FE approximation error. The a posteriori energy norm error estimate and a priori known energy norm errors both converge at optimal rates.

1 Like
4

Thank you both for your quick replies!

@dokken : if a strong solution exists, then the weak solution will coincide with the strong solution. In actual facts, for the numbers I showed I am not even solving any variational problem, I was just interpolating the analytical solution and integrating the strong residual. Sorry if I didn’t make myself clear.

@nate : I am aware [somewhat :smiley:] of a posteriori error estimation, but the thing here is that I am not looking for a way to tell how far my approximate solution is from the real solution, but rather how far the approximate solution to the linearised system is from the approximate solution to the nonlinear problem.
My problem is a(u,v) = (L(u),v), and I solve it as a(u^{n+1},v) = (L(u^n),v). I want to know when to stop by measuring how well a(u^{n+1},v) = (L(u^{n+1}),v) and my initial idea was to do it assembling the L^2 norm of the strong residual.
I wonder if I can use an a posteriori error estimates for this purpose, in a way that gives me an estimate of the error between the numerical solution to the linear problem and the exact solution to the nonlinear problem, but it might be overkill.

After some consideration, I am now thinking I might be able to get away with a discrete solution, that is, after I solve my linear system Au^{n+1} = F(u^n), I can assemble F(u^{n+1}) and compute ||Au^{n+1} - F(u^{n+1})|| and ask that that is below a given tolerance. Does it sound reasonable?

5

What I wrote above applies to the matrix system too. Although posed in a slightly different context, see for example Chapter 2 of Briggs et al. 2000.

6

Got it, thanks a lot for your help, and for the references!

7

Hi @nate,
Sorry for reviving the thread, but I wonder why you only considered the strong residual and not the error related to the flux discontinuity between the elements’ edges. I guess the magnitude of the ‘jump error’ might be smaller than the ‘strong residuals’ error but I’m not sure.

Anyway, if I wish to also account for the flux jumps, following the mentioned paper (Graetsch & Bathe 2005), something like

# `strong' residual
R = div(grad(u)) + f

#Jump discontinuity error
normal = FacetNormal(mesh)
J = jump(grad(u), normal)

c1 = 1.0
c2 = 1.0

a_post_err = sum((c1*cell.circumradius()**2 * assemble_local(R**2*dx, cell) + c2*cell.circumradius() * assemble_local(J**2*dS, cell)) for cell in cells(mesh))**0.5

would make sense in Fenics?

Thanks!

8

My post was so long ago you may choose to assign it to my oversight, laziness or a combination of the two. Definitely place your faith and effort into the work of Graetsch & Bathe 2005 and not my hastily written code.