Clarification on applying Neumann boundary condition or not

8 months ago by
The variational form for a non-linear Poisson's equation

-\nabla \cdot (\epsilon(u) \nabla u) = \rho(u)


F = \int_\Omega{\epsilon(\nabla u \cdot \nabla v) dx} - \int_{\Gamma_N}{\epsilon(\nabla u \cdot \vec{s}) v ds} - \int_\Omega{\rho v dx} = 0

The surface integral applies only to the non-Dirichlet surface boundary because the test function \( v=0\) on Dirichlet boundary. So taking the for example, I tried to modify \(F \) as

F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx
facet_normal = FacetNormal(mesh)
F += -(1+u**2)*inner(grad(u),facet_normal)*v*ds​

in order to explicitly include the boundary integral as given in the weak form. This results in a solution with no Neumann boundary. I tested it and the Neumann surface \(\Gamma_N\) has a non-zero flux across it and that's what I expect.

Adding instead,

F += -(1+u**2)*inner(-grad(u),facet_normal)*v*ds​​

has the same effect as a zero (homogeneous) Neumann boundary.

This essentially imposes the condition \( \nabla u \cdot \vec{s} = -\nabla u \cdot \vec{s} \). Hopefully this will be useful to others as I sometime forget that zero-Neumann does not mean no-Neumann. Is this correct?
Community: FEniCS Project
pde (1) is all messed up...
written 8 months ago by pf4d  
do you mean it does not converge? I changed the script to

# Compute solution
solve(F == 0, u, bcs=bc,
                                          "maximum_iterations": 1000,
written 8 months ago by Chaffra Affouda  
No no, this is fine.  I mean the math of your equation is wrong, there is a dot and nabla in the wrong place.
written 8 months ago by pf4d  
got it. fixed it in the text. thanks.
written 8 months ago by Chaffra Affouda  

1 Answer

8 months ago by
Using the notation from nonlinear Poisson demo, the pde you solve is
- \nabla \cdot \left( \left(1 + u^2 \right) \nabla u \right) = f.
The first example you provide is an example of a bvp with a free-flux boundary condition, meaning that there is no boundary condition imposed along the surface \(\Gamma_N\) so any solution satisfying the pde within \(\Omega\) will suffice.  That is, the flux and value of \(u\) will adjust itself along \(\Gamma_N\) in order to satisfy the variational problem within \(\Omega\).  These boundary conditions are in general not well posed and require care to ensure that the discrete approximations converge.

Regarding your second example, when you change the sign of the surface integral you are not imposing a different boundary condition, but rather solving the same zero-flux-Poisson problem from the demo.  To see why, begin from the variational problem associated with your code; find \(u \subset \mathcal{V} \) which satisfies any Dirichlet boundary conditions on surface \(\Gamma_D\) such that
  \int_{\Omega} \left( 1 + u^2 \right) \nabla u \cdot \nabla v\ d\Omega + \int_{\Gamma_N} \left( 1 + u^2 \right) v \nabla u \cdot n\ d\Gamma_N - \int_{\Omega} f v\ d\Omega = 0 \hspace{10mm} (1)
for all \(v \subset \mathcal{V} \). Using the fact that
  - \int_{\Omega} \nabla \cdot \left( \left( 1 + u^2 \right) \nabla u \right) v\ d\Omega = &+ \int_{\Omega} \left( 1 + u^2 \right) \nabla u \cdot \nabla v\ d\Omega \\ &- \int_{\Gamma_N} \left( 1 + u^2 \right) v \nabla u \cdot n\ d\Gamma_N,
the common weak bilinear volume term can be eliminated resulting in the strong form (in the sense that continuity requirements associated with \(u\) have not been weakend from integration by parts)
  - \int_{\Omega} \nabla \cdot \left( \left( 1 + u^2 \right) \nabla u \right) v\ d\Omega + 2 \int_{\Gamma_N} \left( 1 + u^2 \right) v \nabla u \cdot n\ d\Gamma_N - \int_{\Omega} f v\ d\Omega = 0;
integrating the second-order term of this relation by parts and imposing \(\nabla u \cdot n = 0\) to the surface integral produces
  \int_{\Omega} \left( 1 + u^2 \right) \nabla u \cdot \nabla v\ d\Omega + 2 \int_{\Gamma_N} \left( 1 + u^2 \right) v \nabla u \cdot n\ d\Gamma_N - \int_{\Omega} f v\ d\Omega = 0.
This procedure may be continued indefinitely, and in the limit the integer coefficient of the surface integral goes to infinity.

Therefore, the pde you solve in the second example is the original bvp; Poisson's equation with a zero-flux Neumann boundary.  The extra surface integral does not affect the solution because it is a function of the variable flux which is solved from the linearized system, i.e., it is allowed to vary such that the pde and boundary condition \(\nabla u \cdot n = 0 \) is satisfied uniquely, which in this case implies that the contribution of the extra surface integral --- in the limiting \(k \in \mathbb{N}^+\) case --- is
  \lim_{k \rightarrow \infty} \left\{ k \int_{\Gamma_N} \left( 1 + u^2 \right) v \big[\ \nabla u \cdot n\ \big] \ d\Gamma_N \right\} = \lim_{k \rightarrow \infty} \left\{ k \int_{\Gamma_N} \left( 1 + u^2 \right) v \big[\ 0\ \big] \ d\Gamma_N \right\} = 0.
Further clarifying, you are not implicitly imposing the boundary condition \( \nabla u \cdot n = 0\) by proxy of \( \nabla u \cdot n = - \nabla u \cdot n \), you impose the boundary condition \( \nabla u \cdot n = 0\) and then integrate \( \epsilon \nabla u \cdot n = 0 \) for some \(\epsilon > 0 \) along \(\Gamma_N\).

Here's an example illustrating how the addition of a surface integral proportional to \(\nabla u \cdot n\) along \(\Gamma_N\) to the original Poisson variational problem after integration by parts and imposition of \( \nabla u \cdot n = 0\) does not change the solution (in fact, increasing the magnitude of the coefficient of the surface integral magnifies the numerical error from \(\nabla u^h \cdot n\) with discrete \(u^h \in u \subset \mathcal{V} \)):

from dolfin import *

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
  def inside(self, x, on_boundary):
    return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
File("mesh.pvd") << mesh

V = FunctionSpace(mesh, "CG", 1)

# Define boundary condition
g = Constant(1.0)
bc = DirichletBC(V, g, DirichletBoundary())

# Define variational problem
u  = Function(V)
v  = TestFunction(V)
du = TrialFunction(V)
n  = FacetNormal(mesh)
f  = Expression("x[0]*sin(x[1])", degree=2)
F  = + inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx \
     - Constant(1e10)*inner(grad(u), n)*v*ds

J = derivative(F, u, du)

# Compute solution
solve(F == 0, u, bc, J=J,solver_parameters={"newton_solver":
                                           {"relative_tolerance"  : 1e-6,
                                            "maximum_iterations"  : 200,
                                            "relaxation_parameter": 1.0}})

# Plot solution and solution gradient
plot(u, title="Solution")
plot(grad(u), title="Solution gradient")

# Save solution in VTK format
file = File("nonlinear_poisson.pvd")
file << u

edit 1: provided further mathematical explanations.
edit 2: added further clarifications of math and code.
edit 3: added grammar improvements and more complete description of the code.
edit 4: edited the description of the edit in the last edit to make more sense about that one edit I made way back in the day.

Does this make sense?  
written 8 months ago by pf4d  
I understand your math but that means that dolfin is always trying to satisfy \( \nabla u \cdot n = 0 \) even if the boundary integral is specified in this case. So the boundary integral is adjusted to become essentially 0. I see that's what's happening but I do not quite understand why.
written 8 months ago by Chaffra Affouda  
To help understand you might try looking at the system matrices.
written 8 months ago by pf4d  
Also, I prefer not to think about dolfin doing some math, but rather that dolfin correctly discretizes the math; that is, if you have the math right, dolfin translates your math into code.

This is in contrast to multi-physics-solver packages which actually do the math for you.
written 8 months ago by pf4d  
Your statement that's not clear to me is "Further clarifying, you are not implicitly imposing the boundary condition ∇u⋅n=0 by proxy of ∇u⋅n=−∇u⋅n, you impose the boundary condition ∇u⋅n=0 and then integrate ϵ∇u⋅n=0 for some ϵ>0 along ΓN." Does that mean that dolfin is still applying the  ∇u⋅n=0 boundary condition for the form we are discussing?
written 8 months ago by Chaffra Affouda  
This was to clarify your observation that you were imposing the condition \(\nabla u \cdot n = - \nabla u \cdot n\), which is simply the condition that 1 = -1, implying that \(\nabla \cdot u = 0\).  This is basically noticing that apples are not oranges.  What I have shown is that you could just as easily have determined that \(\nabla u \cdot n = \gamma \nabla \cdot u\) with \(\gamma \neq 0\) for some other choice of added surface integral --- added after you removed the original surface integral resulting from integration by parts --- implying that 1 = \(\gamma\), thus proving that your choice of 1 = -1 was not the reason your pde was applying zero-flux Neumann....  Maybe that is more confusing than it needs to be, but this is the question you asked and that's my best answer.

Regardless, it isn't that dolfin is doing some math, but that you have done some math when you eliminated the surface integral by setting \(\nabla u \cdot n = 0\) by removing the surface integral after your first integration by parts.  Then you integrate over a surface a function which you have already defined as zero.

written 8 months ago by pf4d  
Maybe it would help if you considered how the weak 2nd derivative term without the boundary integral will always implicitly satisfy zero Neumann.  Then, you add some terms to the form which integrate \(\nabla u \cdot n\)...  but there is one and only one choice of surface integral which "fit" the variational problem for non-zero flux; this is the free-flux problem from your first example, where you do not remove the surface integral resulting from integration by parts.
written 8 months ago by pf4d  
I think you could gain insight to why this is so if you look at the proof for the integration-by-parts formula....  also i'm tired and my brain might not be processing things right ...  it seems intuitively that if the surface integral were close to the weak boundary term, it would produce a solution with non-zero flux.  Perhaps you could do a sensitivity report -- vary the coefficient \(\pm\) some amount amount around the true value of the surface integral and plot the coefficient against the norm of the surface flux over \(\Gamma_N\).
written 8 months ago by pf4d  
Ha got it here! That's what was confusing me. I had already subtracted the boundary integral. Thanks for taking the time!
written 8 months ago by Chaffra Affouda  
There is an incredibly rich mathematical theory associated with the Galerkin method; the best source (that I know of) for this theory is the work of Brenner and Scott.
written 8 months ago by pf4d  
[ [ [ [ [ [ bump ] ] ] ] ] ]
written 8 months ago by pf4d  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »