Neumann boundary conditions in NDSolve over nontrivial region

Question

The problem I would like to solve involves diffusion in the following region

reg = ImplicitRegion[-5 <= x <= 5 && -5 <= y <= 5 && x^2 + y^2 >= 1^2, {x, y}];
RegionPlot[reg, AspectRatio -> Automatic]

Additionally, there is an interaction potential $v$ acting in this region:

v = Function[{x, y}, -E^(-x^2 - y^2)];
Plot3D[v[x, y], {x, -5, 5}, {y, -5, 5}, RegionFunction -> Function[{x, y, z}, x^2 + y^2 > 1^2], PlotRange -> {-1, 1}, PlotPoints -> 50]

Physically, I would like to model 2D diffusion confined to this region, in the presence of the biasing potential given by $v$. So what I need to do is solve the PDE $\vec{\nabla} \cdot \left[\vec{\nabla}u(x,y) + u(x,y) \vec{\nabla}v(x,y)\right] = 0$ subject to the boundary conditions

$u(x,-5) = 1$

$u(x,5) =0$

to impose an overall flux, and with zero-flux boundary conditions at the boundary where $x^2 + y^2 = 1$, (as well as on the sides where $x = \pm 5$). So I proceed as follows in Mathematica (10.0.0):

op = Div[Grad[u[x, y], {x, y}] + u[x, y] Grad[v[x, y], {x, y}],{x, y}];
bc = {DirichletCondition[u[x, y] == 1, y == -5],DirichletCondition[u[x, y] == 0, y == 5]};

Note that I did not specify boundary conditions for the circular boundary, nor for the sides $x = \pm 5$, because according to the NDSolve/NeumannValue documentation, these should default to zero-flux boundary conditions.

sol = NDSolve[{op == 0, bc}, u, {x, y} \[Element] reg, Method -> {"FiniteElement","MeshOptions" -> {MaxCellMeasure -> 0.005}}]

The solution looks like this:

Plot3D[Evaluate[u[x, y] /. sol], {x, y} \[Element] reg, PlotRange -> All, AxesLabel -> {"x", "y"}]

Now, this is not quite what I expected, for the following reason:if the PDE is solved correctly, with zero-flux boundary conditions, then the flux into the box (at y = -5) should equal the flux out of the box (at y = 5). You can sort of tell that this is not the case from the plot above, but plotting the normal derivative at the boundaries $y = \pm 5$ makes it clear:

Plot[{Evaluate[(Derivative[0, 1][u][x, 5] + u[x, 5] Derivative[0, 1][v][x, 5]) /.sol],
Evaluate[(Derivative[0, 1][u][x, -5] + u[x, -5] Derivative[0, 1][v][x, -5]) /. sol]}, {x, -5, 5}, PlotRange -> {-0.2, 0}]

Integrating these will give the flux in and the flux out. Clearly the flux in is not equal to the flux out. Since I have checked to make sure that the PDE $\vec{\nabla} \cdot \left[\vec{\nabla}u(x,y) + u(x,y) \vec{\nabla}v(x,y)\right] = 0$ is obeyed in the interior of the region, it seems that I am not correctly setting the Neumann boundary condition at the circular boundary $x^2 + y^2 = 1$. The boundary condition that I want is

$\left[\vec{\nabla}u(x,y) + u(x,y) \vec{\nabla}v(x,y)\right] \cdot \vec{n} = 0$

but I believe that the boundary condition that I'm getting is simply

$\left[\vec{\nabla}u(x,y)\right] \cdot \vec{n} = 0$

How can I impose the correct boundary condition? Have I misunderstood the documentation, which seems to say that my desired boundary condition is the default?

Answer 1

Here is what I think the issue is: Let's look at what NDSolve parses.

Needs["NDSolve`FEM`"]
{state} = 
  NDSolve`ProcessEquations[{op == 0, bc}, u, {x, y} ∈ reg, 
   Method -> {"FiniteElement", 
     "MeshOptions" -> {MaxCellMeasure -> 0.005}}];
femData = state["FiniteElementData"];
femData["PDECoefficientData"]["All"]

{{{{0}}, {{{{0}, {0}}}}}, {{{{{-1, 
      0}, {0, -1}}}}, {{{{0}, {0}}}}, {{{{2 E^(-x^2 - y^2) x, 
      2 E^(-x^2 - y^2) y}}}}, {{4 E^(-x^2 - y^2) - 
     4 E^(-x^2 - y^2) x^2 - 4 E^(-x^2 - y^2) y^2}}}, {{{0}}}, {{{0}}}}

So, there is the diffusion term ({{-1,0},{0,-1}}) the convection term ({2 E^(-x^2 - y^2) x, 2 E^(-x^2 - y^2) y}) and a reaction term. This means NDSolve is modeling something like:

∇⃗ ⋅[∇⃗ u(x,y)] + v(x,y)∇⃗ u(x,y) + u(x,y)∇⃗ v(x,y) = 0

If we manually set the PDE coefficients like this, to be a diffusion coefficient and a conservative convection coefficient (with an added minus sign, as noted in the comments)

methodData = femData["FEMMethodData"];
bcData = femData["BoundaryConditionData"];
sd = state["SolutionData"][[1]];
vd = methodData["VariableData"];

pdeData = 
 InitializePDECoefficients[vd, sd, 
  "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}, 
  "ConservativeConvectionCoefficients" -> {{-{2 E^(-x^2 - y^2) x, 
      2 E^(-x^2 - y^2) y}}}];

This should be more like

∇⃗ ⋅[∇⃗ u(x,y) + u(x,y)∇⃗ v(x,y)]= 0

And discretize and solve the PDE with:

dpde = DiscretizePDE[pdeData, methodData, sd];
dbcs = DiscretizeBoundaryConditions[bcData, methodData, sd];
{l, s, d, m} = dpde["SystemMatrices"];
DeployBoundaryConditions[{l, s, d, m}, dbcs]
lsol = LinearSolve[s, l];
mesh = NDSolve`SolutionDataComponent[sd, "Space"]["ElementMesh"];
sol = {u -> ElementMeshInterpolation[{mesh}, lsol]};

We see that the difference between the in and out flux is relatively small.

Plot[{Evaluate[((Derivative[0, 1][u][x, 5] + 
        u[x, 5] Derivative[0, 1][v][x, 5]) - (Derivative[0, 1][u][
         x, -5] + u[x, -5] Derivative[0, 1][v][x, -5])) /. 
    sol]}, {x, -5, 5}]

Update

In 10.0.2 a better way to input the PDE and get the expected result is made available:

reg = 
  ImplicitRegion[-5 <= x <= 5 && -5 <= y <= 5 && x^2 + y^2 >= 1^2, {x,
     y}];
v = Function[{x, y}, -E^(-x^2 - y^2)];
bc = {DirichletCondition[u[x, y] == 1, y == -5], 
   DirichletCondition[u[x, y] == 0, y == 5]};
α = Grad[v[x, y], {x, y}];
op = Inactive[Div][
   Inactive[Plus][Inactive[Grad][u[x, y], {x, y}], 
    Inactive[Times][α, u[x, y]]], {x, y}];
solI = NDSolve[{op == 0, bc}, u, {x, y} ∈ reg];
Plot[Evaluate[Function[y,
    (Derivative[0, 1][u][x, y] + 
       u[x, y] Derivative[0, 1][v][x, y]) /. 
     solI] /@ {5, -5}], {x, -5, 5}, PlotRange -> {-0.2, 0}]

Answer 2

I would say that you misread the documentation for NeumannValue. I found the documentation fairly difficult to fathom. This was mainly my fault because what they do is substantially different than what I expected. I was trying to cram the round peg that is NeumannValue into the square hole that is my brain.

If you can write the PDE as

Div J = 0

(the form your PDE is in) and you use NeumannValue with its first argument set to zero you are then getting n.J=0. That is if you write your pde as

 Div J = NeumannValue[q - r u,pred]

where u is the dependent variable you are solving for and where pred is true for coordinates that are in the boundaries and false for all other coordinate values then NeumannValue is imposing n.J = q-ru on the boundaries for which pred is true. Thus if you set q =r=0 so that the first argument to NeumannValue is zero then you would in fact be getting the boundary conditions you wanted. If your PDE can't be written as the divergence of a flux (plus time derivatives possibly) or if you want some other boundary conditions then you can't use NeumannValue as far as I know.