Projected gradient different from Neumann BC

[Bratusse]

Hi everyone,
I’m solving Poisson’s equation with a source term defined by an image on a structured mesh. This runs fine (even though there could still be mistakes in my approach), even with a random noise as given in this minimal code. If I’m not mistaken, the way I define the variational problem should impose Neumann conditions on the edge. However by projecting my gradient I get something quite far from du/dn = 0. Is there a way to get a “better” gradient? Here is my minimal code :

from dolfin import *
import numpy as np
from scipy.interpolate import *

h = 60
w = 50
loss = np.random.rand(h,w)
loss_cont = RectBivariateSpline(np.linspace(0.5,h-0.5,h,dtype = np.float64),np.linspace(0.5,w-0.5,w,dtype = np.float64),loss)

# Create mesh and define function space
my_mesh = RectangleMesh(Point(0,0),Point(h,w),h,w)
V = FunctionSpace(my_mesh, 'Lagrange', 1)

vxmap = vertex_to_dof_map(V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Function(V)
f = fem.interpolation.interpolate(f, V)

ks = np.array([i % (w+1) for i in range(V.dim())])
ls = np.array([int(i / (h+1)) for i in range(V.dim())])
f.vector()[vxmap[:]] = loss_cont(ls,ks,grid=False)[:]       # Define source term at its nodes

a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx

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

# Gradient
W = VectorFunctionSpace(my_mesh, "CG", 1)
gradu = project(grad(u), W)
wxmap = vertex_to_dof_map(W)
my_grad = np.array(gradu.vector().get_local()[wxmap]).reshape(h+1,w+1,2)
print(my_grad[:,0])

# Setting the borders
my_grad[:,0,1] = 0
my_grad[:,-1,1] = 0
my_grad[0,:,0] = 0
my_grad[-1,:,0] = 0

# Plot gradient
plot(gradu)

I’m using dolfin on colab, with the package that fem-on-colab provides.

PS : Not sure if I should ask another question in this topic, but I also encountered issues with the display of the solution that, when presented with a few nodes, gives only two colors in the output.

And sometimes it also leaves blank spaces (wasn’t able to reproduce though).

[kamensky]

Note that your problem is not well-posed. A nonzero source term is only compatible with a zero-flux Neumann BC if the integral of the source term is exactly zero. Even with a compatible source term, the solution is only defined up to an arbitrary constant. My answer in the following thread explains these issues in more detail:

White spots in the plot

To expand on @nate’s answer, there are two basic well-posedness issues that come up with pure Neumann BCs. The first is compatibility of the Neumann data with the volume source term. If you have the interior PDE -\nabla\cdot(\nabla u) = f\text{ ,} integrating both sides over \Omega and applying the divergence theorem on the left implies that \int_{\partial\Omega}-\nabla u\cdot\mathbf{n} \,d\partial\Omega = \int_\Omega f\,d\Omega\text{ .} If you also have the Neumann BC of \nabla u\cdot\m…

[Bratusse]

Ah yes you’re right ! Thank you for making me see this error. I guess it could also have an effect on the gradient I’m trying to obtain, I’m gonna try to fix this first and I’ll see !

[Bratusse]

@kamensky Thanks again for the quick reply. I think that my problem’s “well-posedness” has been corrected. I don’t have any more issues with the display. Unfortunately my problem with the gradient persists. It seems to be almost parallel to the edge, but not quite as much as I would like. Here’s the code that I use now :

from dolfin import *
import numpy as np
from scipy.interpolate import *
import matplotlib.pyplot as plt

h = 20
w = 25
loss = np.random.rand(h,w)
loss_cont = RectBivariateSpline(np.linspace(0.5,h-0.5,h,dtype = np.float64),np.linspace(0.5,w-0.5,w,dtype = np.float64),loss)

# Test for PETSc
if not has_linear_algebra_backend("PETSc"):
    info("DOLFIN has not been configured with PETSc. Exiting.")
    exit()

parameters["linear_algebra_backend"] = "PETSc"

# Create mesh and define function space
my_mesh = RectangleMesh(Point(0,0),Point(h,w),h,w)
V = FunctionSpace(my_mesh, 'Lagrange', 1)

vxmap = vertex_to_dof_map(V)

# Define variables
u = TrialFunction(V)
v = TestFunction(V)
f = Function(V)
f = fem.interpolation.interpolate(f, V)

# Define source term at its nodes
ks = np.array([i % (w+1) for i in range(V.dim())])
ls = np.array([int(i / (h+1)) for i in range(V.dim())])
f.vector()[vxmap[:]] = loss_cont(ls,ks,grid=False)[:] - np.mean(loss_cont(ls,ks,grid=False))
print("Sum of source over surface : ", np.sum(f.vector()[:]))

# Define variational problem
a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx

# Assemble system
A = assemble(a)
b = assemble(L)

# Solution Function
u = Function(V)

# Create Krylov solver
solver = PETScKrylovSolver("cg")
solver.set_operator(A)

# Create vector that spans the null space and normalize
null_vec = Vector(u.vector())
V.dofmap().set(null_vec, 1.0)
null_vec *= 1.0/null_vec.norm("l2")

# Create null space basis object and attach to PETSc matrix
null_space = VectorSpaceBasis([null_vec])
as_backend_type(A).set_nullspace(null_space)

null_space.orthogonalize(b);

solver.solve(u.vector(), b)

# Gradient
W = VectorFunctionSpace(my_mesh, "CG", 1)
gradu = project(grad(u), W)
wxmap = vertex_to_dof_map(W)
my_grad = np.array(gradu.vector().get_local()[wxmap]).reshape(h+1,w+1,2)
print("These are supposed to be 0 : ", my_grad[:,0,1])

# Setting the borders
my_grad[:,0,1] = 0
my_grad[:,-1,1] = 0
my_grad[0,:,0] = 0
my_grad[-1,:,0] = 0

# Plot solution and mesh
plot(gradu);

[kamensky]

A few notes:

• Fixing ill-posedness by manipulating the linear algebra problem calls into question the interpretation of the results. I think it’s better practice to formulate the original continuous problem as clearly as possible, rather than trying to fix inconsistencies after discretization. My suggestion (without knowing the full underlying physics) would be to first get compatible data by simply subtracting off the mean value of f from the source term (i.e., f \gets f - \frac{1}{\vert\Omega\vert}\int_\Omega f\,d\Omega), then either use a global Lagrange multiplier to impose some average value for u (demonstrated, e.g., here), or introduce a Dirichlet BC on a single DoF from V, e.g., bc = DirichletBC(V,Constant(0),"x[0]<DOLFIN_EPS && x[1]<DOLFIN_EPS","pointwise"), where "pointwise" is needed to pick off just the corner (since, by default, DOLFIN looks for entire mesh facets matching the geometric criterion).
• In general, when the boundary condition \nabla u\cdot\mathbf{n}=0 is enforced weakly, you would only expect to converge toward it as the mesh is refined. However, the data of your problem also gets rougher with refinement, which may interfere with this. (I’m not sure what numerical analysis exists for introducing random independent noise to the source term at every node; it seems like a questionable practice to me.)
• Your MWE makes some assumptions on the DoF ordering for W that I would not take for granted.
• For some discussion on accurate post-processing of boundary fluxes, you might take a look at the threads here and here.

[Bratusse]

Thanks again for the tips. I will explain my project a little more so that everything becomes clearer. I’m trying to reproduce the results of a paper called Poisson-Based Continuous Surface Generation for Goal-Based Caustics. The general idea is that to compute the surface of the object that will produce the caustics, I have to solve a Poisson equation to which the source term is a loss between the area of the cells of the mesh and the brightness of the pixels of the image that I’m trying to turn into a caustic design. Then, the mesh is updated with the gradient of the solution, so that the areas can converge to the expected brightness, and the process is repeated until the desired precision is achieved. I actually got some pretty good results with the linear algebra solution, but it starts having trouble solving the numerical problem with a large mesh, and probably when some elements start getting flatter.

• I tried the proposed solution of setting a Dirichlet BC, but unfortunately, and as I expected, with the mesh constantly changing it is impossible to keep the integral of f equal to 0 (at least I don’t know how). I am now trying to implement a Lagrange multiplier as in the example, but I’m finding myself having trouble with the mapping of my DoF for the mixed function space, and therefore I can’t define f properly at each node.

• In the code that I provided, the source is simply some random noise so that it can run as is, but in reality it is actually the loss function that I described before. In the end maybe the values of grad(u).n are better the finer the mesh is but I have trouble getting higher than a resolution of 200x200 (trying to get to 512x512), and actually even with a pretty coarse mesh I get some decent results for grad(u) overall.

• I don’t know what MWE is but from the tests that I have conducted it looks like the DoF ordering for W coresponds to how I understood it (and exloited it), but it’s true that I don’t really have proof.

Here is a picture of a simulation that I ran, just for fun

[kamensky]

The paper is really outside my field, but I skimmed it out of curiosity.

• I tried the proposed solution of setting a Dirichlet BC, but unfortunately, and as I expected, with the mesh constantly changing it is impossible to keep the integral of f equal to 0 (at least I don’t know how). I am now trying to implement a Lagrange multiplier as in the example, but I’m finding myself having trouble with the mapping of my DoF for the mixed function space, and therefore I can’t define f properly at each node.

Based on the paper, the constant offset to \phi is unimportant, since you’re only interested in the gradient, so it’s probably easier to just fix \phi at a corner node and correct the mean value of the source term by subtracting a constant from it. The \phi plotted in Figure 4b does look like it has \nabla\phi\cdot\mathbf{n} = 0 on the boundary, which would also be consistent with the edges of the mesh not deforming in Figure 6.

I don’t know what MWE is

It stands for “minimum working example” (as opposed to a whole complicated application code).

I have trouble getting higher than a resolution of 200x200 (trying to get to 512x512), and actually even with a pretty coarse mesh I get some decent results for grad(u) overall.

Speaking of scaling issues, one point that jumped out at me is that the formulation appears to be dimensionally-inconsistent. Every spatial derivative introduces a factor of inverse length to the dimensions, so -\nabla^2\phi = D implies that \phi and D should have different units, but then \nabla\phi is used directly as a replacement for \nabla D. One could argue that it’s a non-physical image processing problem, so everything is dimensionless, but a dimensionally-inconsistent formulation may still lead to unexpected behavior under changes in scale and/or image resolution. Given that \phi is meant to be a smoothed replacement for D, a more typical approach would be to solve a reaction–diffusion (or positive-definite Helmholtz) equation like

-c\nabla^2\phi + \frac{1}{\ell^2}(\phi - D) = 0\text{ ,}

where \ell has units of length and c is a dimensionless constant that controls the strength of the smoothing effect. I’ve seen similar equations used as smoothing filters in other contexts. The additional reaction term also has a nice side effect of eliminating the compatibility requirement between D and Neumann boundary data.

[Bratusse]

Thank you very much for all the time you gave me and for your precious help. My project is at a sufficiently advanced stage for me to move on. I might try your suggestion in the future, it sounds very interesting. In the meantime for anyone who would like to make image-based caustics, here is my tool : https://colab.research.google.com/drive/1OoGreLsagnWwUisVp_rt6sk4dJKr7UdR?usp=sharing