# Laplacian of solution extremely noisy

I have been trying to plot the residual of the solution over my entire domain, but I constantly keep finding that computing the Laplacian of even a well-converged solution ends up creating nonsensically large values. The problem at hand is a system of two Poisson equations with an extra nonlinear forcing term Ku1u2. Here’s the code I used:

import fenics
import numpy as np

ndims = 2
domain_extent = [[-2.0,2.0],[-2.0,2.0]]
p_bottomleft = fenics.Point(*[extent[0] for extent in domain_extent])
p_topright = fenics.Point(*[extent[1] for extent in domain_extent])
mesh = fenics.RectangleMesh.create(p=[p_bottomleft,p_topright],n=[100 for _ in range(ndims)],cell_type=fenics.CellType.Type.triangle)

P1 = fenics.FiniteElement('CG', fenics.triangle, 2)
element = fenics.MixedElement([P1 for _ in range(ndims)])
FuncSpace = fenics.FunctionSpace(mesh, element)
TestFuncs = fenics.TestFunctions(FuncSpace)
SolnFuncs = fenics.Function(FuncSpace)
u = fenics.split(SolnFuncs)

tol = 1e-14
bc1l = fenics.DirichletBC(FuncSpace.sub(0), fenics.Constant(-0.15), lambda x, ob: ob and fenics.near(x[0],-2.0, tol))
bc1r = fenics.DirichletBC(FuncSpace.sub(0), fenics.Constant(0.15), lambda x, ob: ob and fenics.near(x[0],2.0, tol))
bc1t = fenics.DirichletBC(FuncSpace.sub(0), fenics.Constant(0.30), lambda x, ob: ob and fenics.near(x[1],-2.0, tol))
bc1b = fenics.DirichletBC(FuncSpace.sub(0), fenics.Constant(0.30), lambda x, ob: ob and fenics.near(x[1],2.0, tol))
bc2l = fenics.DirichletBC(FuncSpace.sub(1), fenics.Constant(0.30), lambda x, ob: ob and fenics.near(x[0],-2.0, tol))
bc2r = fenics.DirichletBC(FuncSpace.sub(1), fenics.Constant(0.30), lambda x, ob: ob and fenics.near(x[0],2.0, tol))
bc2t = fenics.DirichletBC(FuncSpace.sub(1), fenics.Constant(-0.15), lambda x, ob: ob and fenics.near(x[1],-2.0, tol))
bc2b = fenics.DirichletBC(FuncSpace.sub(1), fenics.Constant(0.15), lambda x, ob: ob and fenics.near(x[1],2.0, tol))
bcs = [bc1l,bc1r,bc1t,bc1b,bc2l,bc2r,bc2t,bc2b]

diffusion_coeff = 10.0
source_coeff = 1.0
forcing_expression_1 = fenics.Expression('0.33*sin(3.1415*(x[0]-0.8*x[1])/4.0)', degree = 3)
forcing_expression_2 = fenics.Expression('-0.175*cos(3.1415*x[0]*x[1]/4.0)', degree = 3)

eq1_linear_part = \
- forcing_expression_1 * TestFuncs[0] * fenics.dx
eq1 = eq1_linear_part + source_coeff * u[0] * u[1] * TestFuncs[0] * fenics.dx
eq2_linear_part = \
- forcing_expression_2 * TestFuncs[1] * fenics.dx
eq2 = eq2_linear_part + source_coeff * u[0] * u[1] * TestFuncs[1] * fenics.dx

solver_parameters = {
"newton_solver":{"maximum_iterations":50, 'relative_tolerance': 1e-6, 'absolute_tolerance': 1e-5, "error_on_nonconvergence": False, "linear_solver": "lu"}
}

SolnFuncs.vector().set_local(2*np.random.rand(SolnFuncs.vector().size())-1.0)
SolnFuncs.vector().apply("")

fenics.solve(eq1 + eq2 == 0, SolnFuncs, bcs, solver_parameters = solver_parameters)
#plot the solution
import matplotlib.pyplot as plt
p = fenics.plot(SolnFuncs[0])
plt.colorbar(p)
plt.savefig('/build/tmp/u0.png')
plt.figure()
#compute laplacian of u0 and plot it
p = fenics.plot(lFF0)
plt.colorbar(p)
plt.savefig('/build/tmp/u0_laplacian.png')
plt.figure()


Below are the solution and the corresponding Laplacian (see 1st reply to this post for the Laplacian):

As you can see below, the Laplacian has very large spikes at the corners and the values in the rest of the domain, I have found, do not really lead to mean residual values anywhere near the value printed to stdout during the solution process. Is there a way to improve the quality of the values computed?

p.s. this was computed with fenics installed from the package available via apt on Ubuntu 20.04

The forums didn’t allow me to upload 2 pictures so here’s the plot of the Laplacian of u0:

Note that the Laplacian of a CG 2 function is not a CG 2 function.
The gradient of a CG 2 function is in the vector DG 1 space.
The trace of this function is constant inside each cell, but not well defined between cells.

1 Like

Which element would you suggest for this type of application, in that case? I noticed that higher order CG elements result in even worse behaviour than CG 2. Naturally CG 1 has 0 laplacians everywhere.

Alternatively, is there a better/recommended way to get the cell-wise residuals?

I guess I would use a DG 0 space, but I am not sure why you want to visualize this. The reason for using finite elements is that you reduce continuity of the solution space.

I need the cell-wise residuals to integrate it into a reinforcement learning style control scheme as an objective for the model.

I get the following error if I use the DG 0 space for the solution variable:

ufl.log.UFLException: This integral is missing an integration domain.

Is there a way to e.g. solve the equation for u1 and u2 using CG 2 elements but then compute the laplacian using some other type of element after projecting the solutions to another function space? As a very naive solution I tried the following:

P2 = fenics.FiniteElement('DG', fenics.triangle, 0)
FuncSpace_Laplacian = fenics.FunctionSpace(mesh, element)
lFF0 = fenics.project(u[0], FuncSpace_Laplacian)
p = fenics.plot(lFF0_values)


which did not work unfortunately, with the following error:

Traceback (most recent call last):
File “dataset/residual_question.py”, line 61, in
lFF0 = fenics.project(u[0], FuncSpace_Laplacian)
File “/usr/lib/python3/dist-packages/dolfin/fem/projection.py”, line 129, in project
L = ufl.inner(w, v) * dx
File “/usr/lib/python3/dist-packages/ufl/operators.py”, line 169, in inner
return Inner(a, b)
File “/usr/lib/python3/dist-packages/ufl/tensoralgebra.py”, line 161, in new
error(“Shapes do not match: %s and %s.” % (ufl_err_str(a), ufl_err_str(b)))
File “/usr/lib/python3/dist-packages/ufl/log.py”, line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shapes do not match: and .

I think you’re missing a fundamental point that @dokken is making. The broken Laplace operator (\nabla^2 defined on each element)

\nabla^2_h : {V}^{h,2} \rightarrow {V}_\text{DG}^{h,0}.

Your MWE is rather confusing so I’ll just write out something simpler

from dolfin import *

mesh = UnitSquareMesh(32, 32)
x = SpatialCoordinate(mesh)

V = FunctionSpace(mesh, "CG", 2)
VDG0 = FunctionSpace(mesh, "DG", 0)

u = TrialFunction(V)
v = TestFunction(V)

L = sin(pi*x[0])*sin(pi*x[1])*v*dx

bcs = DirichletBC(V, Constant(0.0), "on_boundary")

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

import matplotlib.pyplot as plt
plt.title("$u_h(x,y)$")
plt.title(r"$\Pi_{V_{DG}^{0,h}} (\nabla^2 u_h(x,y))$")