How to calculate the Laplacian of a known solution correctly?

Legacy DOLFIN
1

Hello. I am trying to make a function that takes a known solution and returns its Laplacian. The sig expression is there to make my 1D solution match the 3D one which is spherically symmetric. for a specified coordinate system. When I run the code I get the correct solution (Laplacian=6) but with large oscillations at the boundaries (see figure). I have also confirmed this is also a problem when using a 2D mesh also.

Am I using the a wrong approach to calculating the Laplacian? I know I could use something like ‘project(div(grad(f)))’ but that would only work in Cartesian coordinates, where as using the integration by parts allows me to works with different sig.

FEinCS version - 2019.1.0 (installed with conda-forge).

Thanks in advance.

import dolfin as d

mesh = d.UnitIntervalMesh(1000)
sig = d.Expression("pow(x[0], 2)", degree=0)  # in spherica coordinates.

V = d.FunctionSpace(mesh, 'CG', 1)
v = d.TestFunction(V)
u = d.TrialFunction(V)

f = d.interpolate(d.Expression("pow(x[0], 2)", degree=1), V)

solver = d.KrylovSolver(method="richardson", preconditioner="icc")

n = d.FacetNormal(mesh)
ds = d.Measure('ds', domain=mesh)

A = d.assemble(u*v*sig*d.dx)
G = d.assemble(d.dot(d.grad(f), n)*v*sig*ds
               - d.dot(d.grad(f), d.grad(v*sig))*d.dx)

laplacian = d.Function(V)
solver.solve(A, laplacian.vector(), G)

d.plot(laplacian)

Laplacian_calc
Laplacian_calc370×248 3.75 KB

2

It looks like you’re solving
\sigma u = \sigma\Delta f
Are you sure this is what you need? What is then still the point of \sigma?

It should be Laplacian=2 though?

It seems like there is a bunch of small inconsistencies at play

implies that it is a piecewise constant. Since it is inside a gradient, this is harmful. You’ll need

sig = d.Expression("pow(x[0], 2)", degree=2) 
sig = d.interpolate(sig, V)

But you’re now interpolating a quadratic function on a linear function space, so that’ll create some noise on your result.

Then,

f = d.interpolate(d.Expression("pow(x[0], 2)", degree=1), V)

should equally be:

f = d.interpolate(d.Expression("pow(x[0], 2)", degree=2), V)

Furthermore, your interative solver is throwing in some round off errors, so I’d simple use d.solve.

Then, the evaluation of \nabla f \cdot n on the boundary is a ‘variational crime’, which has a significant impact due to the interpolation of f onto a low-order space (linears). Would this be possible in your usecase?

df = d.interpolate(d.Expression("2*x[0]", degree=1), V)
G = d.assemble(df*n[0]*v*sig*ds
               - d.dot(d.grad(f), d.grad(sig*v))*dx)

But then we’re still not there… as there is still a missmatch on the left side. I assume this is due to sig being zero there, and the solution being ill defined. Not sure if there is a proper way around. (Maybe this is not what you actually mean to solve? See my first sentence.)

So that would result in:

import dolfin as d
from matplotlib import pyplot as plt

mesh = d.UnitIntervalMesh(100)
sig = d.Expression("pow(x[0], 2)", degree=2)  # in spherica coordinates.

V = d.FunctionSpace(mesh, 'CG', 1)
v = d.TestFunction(V)
u = d.TrialFunction(V)

f = d.interpolate(d.Expression("pow(x[0], 2)", degree=2), V)
sig = d.interpolate(sig, V)

n = d.FacetNormal(mesh)
ds = d.Measure('ds', domain=mesh)
dx = d.Measure('dx', domain=mesh)

A = d.assemble(u*v*sig*dx)

df = d.Expression("2*x[0]", degree=1)
G = d.assemble(df*n[0]*v*sig*ds
               - d.dot(d.grad(f), d.grad(sig*v))*dx)

laplacian = d.Function(V)
d.solve(A, laplacian.vector(), G)

d.plot(laplacian)
plt.show()

Lastly (and probably least usefully), you would get the near exact result for

V = d.FunctionSpace(mesh, 'CG', 2) # Second order
2 Likes
3

Thank you for your response.

As I mentioned in my post, the problem is spherically symmetric. Hence I need to defined the problem in spherical coordinates. \sigma = r^2 comes from the volume element dV = r^2 sin^2(\theta) dr d\theta d\phi \rightarrow 4 \pi \sigma dr, where in the code I divided out the common factor of 4 \pi from each integral. This is required so the solution on my 1D mesh matches a 3D result with this imposed symmetry.

On this note since we are solving the Laplacian in spherical coordinates the answer is 6 :

\nabla^2 (r^2) = \partial_r^2 r^2 + \frac{2}{r} \partial r^2 = 2 + 4 = 6.

Apologies, sig was originally outside the grad (i.e. d.dot(d.grad(f), d.grad(v))*(sig)*dx). Are there any problems with this?

That said I know the problem with the boundary values are unrelated to sig as in the case of using Cartesian coordinates (sig = 1) I still have the same problem.

Running your code I got the wrong answer, but moving the sig out of the grad and changing the degree of V to 2, as you suggested, did give me the right answer with an error of 1e-10. I also found lowering the degree of sig had no effect (is this ill advised?). I did find a problem though I wanted to check. If I use a very long mesh the error accumulates for large x[0] and increasing the mesh resolution also increases this error. Is this something that can be avoided?

Code:

import numpy as np
import dolfin as d
from matplotlib import pyplot as plt

deg_V = 2

# mesh = d.UnitIntervalMesh(100)
mesh = d.RectangleMesh(d.Point(0.0, 0.0), d.Point(1000.0, 1.0), 10000, 10)

V = d.FunctionSpace(mesh, 'CG', deg_V)
v = d.TestFunction(V)
u = d.TrialFunction(V)

f = d.interpolate(d.Expression("pow(x[0], 2)", degree=deg_V), V)
# sig = d.interpolate(d.Expression("pow(x[0], 2)", degree=0), V)
sig = d.Constant(1) # changed to Cartisan coordinates.

n = d.FacetNormal(mesh)
ds = d.Measure('ds', domain=mesh)
dx = d.Measure('dx', domain=mesh)

A = d.assemble(u*v*sig*dx)
G = d.assemble(d.dot(d.grad(f), n)*v*sig*ds
               - d.dot(d.grad(f), d.grad(v))*sig*dx)

laplacian = d.Function(V)
d.solve(A, laplacian.vector(), G)

# d.plot(laplacian)

x = np.linspace(0, 1000, 10000)
Lp = [laplacian(xi, 0) for xi in x]

Laplacian_calc(2)
Laplacian_calc(2)384×259 6.99 KB

4

That looks more safe. Just to make sure you’re doing what you intend to be doing; when you move \sigma our of the gradient, you’re effectively solving a different problem:
\sigma u = \nabla\cdot(\sigma \nabla f)

Regarding the drift in your results, I actually don’t know… Given the simplicity and exactness of the toy problem, I’d have expected you should get much closer to machine accuracy. But I can’t spot any mistakes.

5

Yes that is the correct equation since it reproduces the Laplacian in my different coordinates.

As for the drift, I have tried various things (changing the mesh, changing the initial value of laplacian, changing the degree) and nothing has worked.

I don’t know if this might help determine where the proplem is but I found by testing different functions (i.e. f(x)=(x, x^2, x^3)) that the error scales with f(x) (i.e. \epsilon = (\nabla^2 f - ans)/f = 10^{-11} + noise). I have an example plot for the f(x)=x^2 case below.

rel_error
rel_error399×262 5.89 KB

I suspect it might come from d.grad(f), since if I replace G with G = d.assemble(d.grad(f)[0]vdx) I get an equivalent scaling relationship (see below code).

import numpy as np
import dolfin as d
from matplotlib import pyplot as plt

deg_V = 2

mesh = d.RectangleMesh(d.Point(0.0, 0.0), d.Point(100.0, 1.0), 1000, 10)

V = d.FunctionSpace(mesh, 'CG', deg_V)
v = d.TestFunction(V)
u = d.TrialFunction(V)

f = d.interpolate(d.Expression("pow(x[0], 2)", degree=deg_V), V)

n = d.FacetNormal(mesh)
ds = d.Measure('ds', domain=mesh)
dx = d.Measure('dx', domain=mesh)

A = d.assemble(u*v*dx)
G = d.assemble(d.grad(f)[0]*v*dx)  # Just calculating grad.

y = d.Function(V)
d.solve(A, y.vector(), G)

x = np.linspace(1e-3, 100, 10000)
ep = [abs(y(xi, 0)-(2*xi))/f(xi, 0.5) for xi in x]

plt.figure()
plt.ylabel(r"$\epsilon$")
plt.xlabel('x')
plt.yscale('log')
plt.plot(x, ep)

rel_error_grad
rel_error_grad399×262 7.94 KB

1 Like
6

After further checking it is definitely a problem coming from the ‘grad’ function. (the error is unaffected by changing how I assign the values of f, changing the degree, or by calculating the grad using a vector space). As such I don’t know how to address it or whether this changes the nature of the topic?

import numpy as np
import dolfin as d
from matplotlib import pyplot as plt

mesh = d.RectangleMesh(d.Point(0.0, 0.0), d.Point(100.0, 1.0), 1000, 10)

V = d.FunctionSpace(mesh, 'CG', 2)
v = d.TestFunction(V)
u = d.TrialFunction(V)

f = d.interpolate(d.Expression("x[0]", degree=0), V)
y = d.project(d.grad(f)[0], V)

x = np.linspace(1e-3, 100, 10000)
Dy = [abs(y(xi, 0)-1) for xi in x]
plt.plot(x, Dy)

error_grad
error_grad362×259 6.79 KB