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
Wthat 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.