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 K*u1*u2. 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 = \
+ fenics.inner(fenics.grad(u[0]),fenics.grad(TestFuncs[0]))*diffusion_coeff*fenics.dx \
- forcing_expression_1 * TestFuncs[0] * fenics.dx
eq1 = eq1_linear_part + source_coeff * u[0] * u[1] * TestFuncs[0] * fenics.dx
eq2_linear_part = \
+ fenics.inner(fenics.grad(u[1]),fenics.grad(TestFuncs[1]))*diffusion_coeff*fenics.dx \
- 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
lFF0 = fenics.project(fenics.div(fenics.grad(u[0])), FuncSpace.sub(0).collapse())
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