Neumann boundary condition strange behaviour

Question

Hi all,

I am trying to solve a PDE with Neumann boundary conditions on the unit cube. I can provide the PDE itself if necessary, but I think the problem already becomes clear from the variational formulation which I used. Here I consider the first component as time and the other two as space. I impose zero Neumann boundary condition in space and a some function in time. This seems to run as I expect and the solution is fine. But if I then compute the gradient of this I get very strange results and the initial condition is not satisfied on the boundary anymore. It looks like it is blurred out and slightly shifted which really surprises me. I used Paraview to visualise the results and to compare them. Is anyone able to comment on that?

Here is my code:

Loading all relevant packages.

from fenics import *

Relevant constants and initial functions

delta = Constant(0.001)
g1 = Expression('exp(-150pow(x[1]-0.4,2) - 150pow(x[2]-0.4,2))', degree=0)
g2 = Expression('exp(-150pow(x[1]-0.6,2) - 150pow(x[2]-0.6,2))', degree=0)

Create Boundary conditions for PDE

Boundary for t=0

class Boundary0(SubDomain):

def inside(self, x, on_boundary):
# return True if t=0
    return bool(near(x[0], 0) and  on_boundary)

Boundary for t=1

class Boundary1(SubDomain):

def inside(self, x, on_boundary):
# return True if t=1
    return bool(near(x[0], 1) and  on_boundary)

Initializing the mesh,ds and functions

mesh = UnitCubeMesh(12,20,20)
V = FunctionSpace(mesh, "CG", 1)
V0 = FunctionSpace(mesh, "DG", 0)
W = VectorFunctionSpace(mesh, "DG", 0 , 3)#,constrained_domain=PeriodicBoundary()
mu = interpolate(Expression(('0', "0", "0"), degree=0), W)
g1 = interpolate(g1,V0)
g2 = interpolate(g2,V0)
boundary_markers = FacetFunction('size_t', mesh,0)
NB0=Boundary0()
NB0.mark(boundary_markers, 1)
NB1=Boundary1()
NB1.mark(boundary_markers, 2)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

PDE which we want so solve in variational form

v = TestFunction(V)
phi = TrialFunction(V)
a = (delta/2phiv+ dot(grad(phi)/2., grad(v)))dx
b = -(g1v)ds(1) + (g2v)ds(2) - dot(mu, grad(v))dx
phi = Function(V)
solve(a==b,phi)

Checking boundary condition

mu = project(grad(phi),W)

Save files

musol, m1, m2 = mu.split(deepcopy=True)
vtkfile = File('Test/solution.pvd')
vtkfile << musol
vtkfile_1 = File('Test/initial1.pvd')
vtkfile_1 << g1
vtkfile_2 = File('Test/initial2.pvd')
vtkfile_2 << g2

It would be great if someone can explain this behaviour. I am really confused and can not explain the effect at all.

Thanks,
Sebastian

Answer 1

Hi there,

the line

phi = Function(V)

seems odd, because you have already defined phi above as a trial function. Why do you need that line?

Quoting from the FEniCS tutorial, I would attempt to solve your problem along these lines (making the obvious adjustments):

    # Define variational problem
     u = TrialFunction(V)
     v = TestFunction(V)
     f = Constant(-6.0)
     a = dot(grad(u), grad(v))*dx
     L = f*v*dx

    # Compute solution
    u = Function(V)
    solve(a == L, u, bc)

Does that help?

Comments

Isn't this example doing exactly the same? First define the variational problem with trial functions and then setting it to a function after?

---

Yes, I think you are right.

How sure are you about the CG/ DG setup of your function spaces? Is that standard? Also, you project phi, an element of the CG space V, to the DG space W, which might be non-trivial or not even well defined.

---

The idea of this setup is the following. If I compute the derivative of a piecewise linear function I would expect it to be piecewise constant. However, the gradient of a fenics function is not a regular function and you have to project it before you can use it as one again. At least this is my experience with it so far. With the above reasoning, I would expect that the function stays the same under projection, only the type of it changes to a fenics function. Does this make any sense?

---

OK, but why do you need a DG function space? Also, why do you assume that your solution phi will be a piecewise linear function?

---

Because my space of elements in which I solve the problem is 'CG' 1. So at least my numerical solution phi is continuous, piece wise linear then. And after computing the derivative it will be piece wise constant then. So the gradient of phi should be in 'DG' 0 I would assume...

---

Sorry for my late reply!

Ok, that also makes sense.

I currently can't spot any more issues with your code. A general comment would be to make sure your problem is a well-defined PDE problem in the mathematical sense (see literature) to make sure a solution with the required regularity actually exists.

For example, you have a small constant delta in your code that enters the definition of your linear form a(phi, v). Maybe increase that constant to see whether the problem persists, because you might have a boundary layer somewhere which are hard to deal with numerically.