How to check that the Neumann boundary condition is correctly applied

variational formulation
1

Hi everyone,
I want to check (numerically) whether or not the Neumann boundary condition is correctly applied (after I get a solution) for the following problem.
question fenics

Here, Ix, Iy, It are image derivatives and u, v are velocity components.
Below is my main code.

(nx, ny) = I.shape
mesh = RectangleMesh(Point(0, 0), Point(ny-1, nx-1), ny-1, nx-1)
x = mesh.coordinates()
V = FunctionSpace(mesh, "CG",1)
f_x = image_to_CG1(Ix, V)# CG 1 interpolation
f_y = image_to_CG1(Iy, V)
f_t = image_to_CG1(It, V)
W = VectorFunctionSpace(mesh, "CG", 1, 2)

w = TrialFunction(W)
q = TestFunction(W)
(u, v) = split(w)
(q1, q2) = split(q)

lamda = Constant(0.01)

a = lamda*dot(grad(u), grad(q1))*dx + \
    lamda*dot(grad(v), grad(q2))*dx + \
    (f_x*u+f_y*v)*f_x*q1*dx + (f_x*u+f_y*v)*f_y*q2*dx

L = - f_t*f_x*q1*dx - f_t*f_y*q2*dx

c = Constant((1.0, 1.0))

w = Function(W)
init = interpolate(c, W) #nonzero_initial_guess
w.assign(init)
A = assemble(a)
b = assemble(L)
max_iter = 10000
solver = KrylovSolver('gmres', 'ilu')#cg/gmres  jacobi/ilu
solver.parameters['absolute_tolerance'] = 1e-7
solver.parameters['relative_tolerance'] = 1e-5
solver.parameters['maximum_iterations'] = max_iter
#print(solver.parameters.keys())

solver.parameters['nonzero_initial_guess'] = True
solver.solve(A, w.vector(), b)

Thank you for any feed back.

2

I think you should be able to calculate the flux on \delta\Omega via:

ds = Measure('ds', domain=mesh)
n = FacetNormal(mesh)
flux_u = assemble( dot(n, u)*ds)
print(flux_u)

Let me know if this works for you. Consider that it might not be exactly 0 due to the approximation with your mesh, however it should slowly converge to zero for finer meshes.

3

I have just used your code and I got the value: 556.587629374805 :unamused:
Any comments.

4

Consider the following example:

from dolfin import *

mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 2)
u, v = TrialFunction(V), TestFunction(V)

a = dot(grad(u), grad(v))*dx
L = Constant(1.0)*v*dx

LEFT, RIGHT, TOP, BOTTOM = 1, 2, 3, 4
ff = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
AutoSubDomain(lambda x: x[0] < DOLFIN_EPS).mark(ff, LEFT)
AutoSubDomain(lambda x: x[0] > 1.0 - DOLFIN_EPS).mark(ff, RIGHT)
AutoSubDomain(lambda x: x[1] < DOLFIN_EPS).mark(ff, BOTTOM)
AutoSubDomain(lambda x: x[1] > 1.0 - DOLFIN_EPS).mark(ff, TOP)

bcs = [DirichletBC(V, Constant(0.0), ff, BOTTOM)]

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

n = FacetNormal(mesh)
ds = Measure("ds", subdomain_data=ff)
solution_trace_norm = assemble(dot(grad(uh), n)**2*ds((LEFT, RIGHT, TOP)))**0.5
print(f"Homogeneous Neumann BC trace: {solution_trace_norm}")

which gives me:

Homogeneous Neumann BC trace: 1.0071881730182207e-13
5

@nate Would the following be correct for how the author has written their code?

ds = Measure('ds', domain=mesh)
n = FacetNormal(mesh)
trace_x = assemble(dot(n, grad(w.sub(0)))**2*ds)**0.5
trace_y = assemble(dot(n, grad(w.sub(1)))**2*ds)**0.5
6

@Sven It gives me
7.11056026957267
5.808202663765562

7

I’m not sure what your image_to_CG1 is doing but consider the following MWE. Perhaps tightening your tolerance on the KrylovSolver will help based on your coefficients.

from fenics import *

(nx, ny) = (120, 120)
mesh = RectangleMesh(Point(0, 0), Point(ny-1, nx-1), ny-1, nx-1)
x = mesh.coordinates()
V = FunctionSpace(mesh, "CG",1)
f_x = Constant(1)# CG 1 interpolation
f_y = Constant(1)
f_t = Constant(1)

W = VectorFunctionSpace(mesh, "CG", 1, 2)
w = TrialFunction(W)
q = TestFunction(W)
(u, v) = split(w)
(q1, q2) = split(q)

lamda = Constant(0.01)

a = lamda*dot(grad(u), grad(q1))*dx + \
    lamda*dot(grad(v), grad(q2))*dx + \
    (f_x*u+f_y*v)*f_x*q1*dx + (f_x*u+f_y*v)*f_y*q2*dx

L = - f_t*f_x*q1*dx - f_t*f_y*q2*dx

c = Constant((1.0, 1.0))

w = Function(W)
init = interpolate(c, W) #nonzero_initial_guess
w.assign(init)
A = assemble(a)
b = assemble(L)
max_iter = 10000
solver = KrylovSolver('gmres', 'ilu')#cg/gmres  jacobi/ilu
solver.parameters['absolute_tolerance'] = 1e-8
solver.parameters['relative_tolerance'] = 1e-6
solver.parameters['maximum_iterations'] = max_iter
solver.parameters['nonzero_initial_guess'] = True
solver.solve(A, w.vector(), b)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0, 1e-14)

edge_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
edge_markers.set_all(0)
left = Left()
left.mark(edge_markers, 1)

V = FunctionSpace(mesh, 'CG', 1)
ds = Measure('ds', domain=mesh, subdomain_data=edge_markers)
n = FacetNormal(mesh)
print(assemble(dot(n, grad(w.sub(0)))**2*ds)**0.5)
print(assemble(dot(n, grad(w.sub(1)))**2*ds)**0.5)

The above yields the following:

5.269415601331755e-06
5.269415617215526e-06

Let me know if restricting the solver parameters helps.

8

@nate This above code yields:
Homogeneous Neumann BC trace: 9.18124640398017

9

@Sven image_to_CG1(Ix, V) does dof ordering of the nx*ny matrix Ix and then interpolate it into V to get a function f_x.

11

@nate Should I expect to get zero flux on the boundary even if the solution is not quite good?