Laplacian of solution extremely noisy

mathematics
1

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 = \
    + 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):

u0

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

2

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

3

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
4

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?

5

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.

6

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

7

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)
lFF0_values = fenics.project(fenics.div(fenics.grad(lFF0)), 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 .

8

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)

a = dot(grad(u), grad(v))*dx
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)

laplace_u = project(div(grad(u)), VDG0)

import matplotlib.pyplot as plt
plt.gcf().add_subplot(1, 2, 1)
plot(u)
plt.title("$u_h(x,y)$")
plt.gcf().add_subplot(1, 2, 2)
plot(laplace_u)
plt.title(r"$\Pi_{V_{DG}^{0,h}} (\nabla^2 u_h(x,y))$")
plt.show()

image

1 Like
9

Thank you for the instructive example. For my own problem, I was able to fix the issue but I had to use the CG3 space for the solution variables and then project to DG0 for the Laplacian - using CG2 and DG0 still ended up being problematic. I am a bit unsure why the added nonlinear term causes that - I suppose I need to brush up on my FEM knowledge a little bit to get to the bottom of the issue.