NeumannValue is producing incorrect results

Question

Cross-posted in Wolfram community, an online Mathematica notebook is available there.

---

The problem under investigation is solving the Laplacian equation to model seepage. See the image below for the region and the model. The function u[x,y] represents the total head at the coordinate [x,y]. The region under study is a rectangle. There are two Dirichlet boundary conditions and three Neumann boundary conditions in the model. In the Neumann boundary conditions, there is no flux through the boundary edges on the left, the right and the bottom. That is to say that the partial derivative of u[x,y] with respect to y is zero when y==0. In the meantime, the partial derivative of u[x,y] with respect to x is zero when x==0 or x is the right edge of the region. The screenshot shows the correct solution in a FEM program.

However, I defined the problem in mathematica and obtained incorrect results. I checked the gradient vector of the solution.

v = -Grad[solution[x, y], {x, y}]

This is the result of v. It is clear that there is "FLUX" through the Neumann boundaries. There should be no flux through the identified domain.

How can I fix this error?

The Mathematica notebook is can be found herein. Notebook

Here's the complete code sample:

eqn = Laplacian[u[x, y], {x, y}] == 0;
width = 48;
height = 12;
Ω = Rectangle[{0, 0}, {width, height}];
Subscript[Ω, 1] = Polygon[{{0, 0}, {43, 0}, {23, 10}, {20, 10}}];
dcond1 = DirichletCondition[u[x, height] == 22, x <= width/3];
dcond2 = DirichletCondition[u[x, height] == 18, x > (2 width)/3];

solution = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 
  NeumannValue[0, y == 0] + NeumannValue[0, x == 0] + 
  NeumannValue[0, x == width] + 
  NeumannValue[0, y == height && (width/3 < x < (2 width)/3)], {dcond1, dcond2}},
    u, {x, y} ∈ Ω];

ContourPlot[solution[x, y], {x, y} ∈ Ω,
 AspectRatio -> 1/3, PlotRange -> All, AxesLabel -> {"x", "y"}]

v = Grad[-solution[x, y], {x, y}];
VectorPlot[v, {x, y} ∈ Ω, AspectRatio -> 1/3, 
 PlotRange -> All, AxesLabel -> {"x", "y"}, PlotLegends -> Automatic]
Plot3D[solution[x, y], {x, y} ∈ Ω,
 AspectRatio -> 2/3, PlotRange -> All, AxesLabel -> {"x", "y", "u(x,y)"}, 
 PlotTheme -> "Business"]

Answer 1

This is a bit funny. It's not NeumannValue that's causing trouble. It's all because you've specifying the DirichletCondition in a mixed manner. You specify constraint of x in 2nd argument, but include constraint of y in 1st argument as u[x, height]. Though not clearly stated in the document, this syntax is not valid (at least for now). If you check the examples in document of DirichletCondition, you'll find constraints are always set in 2nd argument.

"But Mathematica doesn't complain! This suggests I'm not setting DirichletCondition in wrong way!" This is the funny part. The error checking of DirichletCondition turns out to be rather rough. Even with

dcond1 = DirichletCondition[u[arguments, aren' t] == 22, x <= width/3];
dcond2 = DirichletCondition[u[checked, at all…] == 18, x > (2 width)/3];

NDSolve silently produces the incorrect result as shown in the question.

So the solution is simple: set the Dirichlet conditions in documented way.

dcond1 = DirichletCondition[u[x, y] == 22, y == height && x <= width/3];
dcond2 = DirichletCondition[u[x, y] == 18, y == height && x > (2 width)/3];

Now the solution is desired:

eps = 1;
ContourPlot[solution[x, y], {x, y} ∈ Ω, AspectRatio -> 1/3, 
  ColorFunction -> "TemperatureMap"]~Show~
 VectorPlot[v, {x, 0, width}, {y, eps, -eps + height}(*{x,
  y}∈Ω*), AspectRatio -> 1/3
  ]