Incorrect solution of 3D PDE with Neumann conditions?

Question

I would like to solve the following PDE: $$ \Delta \mu (x,y,z)=\mu$$ Over the parallelepiped: $x \in [0,100] $, $y \in [0,100] $ and $z \in [-10,0] $

With the following Neumann boundary conditions: $$ \mathbf n \cdot \nabla \mu |_{z=0}= x y, \ \ \ \ \ \ \ \mathbf n \cdot \nabla \mu |_{ rest \ of \ boundary}= 0$$

---

I use Mathematica 11.0. First, I specify the region:

reg = Parallelepiped[{0, 0, -10}, {{100, 0, 0}, {0, 100, 0}, {0, 0, 10}}]

Then I solve the diffusion equation:

sol = NDSolveValue[Laplacian[mu[x, y, z], {x, y, z}] - mu[x, y, z] == NeumannValue[- x y, z == 0]},
    mu, {x, y, z} \[Element] reg]

I get a solution (however, I get the mistake NDSolveValue::femibcnd: No DirichletCondition or Robin-type NeumannValue was specified for {mu}; the result may not be unique.). Afterwards, I want to check whether the second boundary condition is satisfied (in documentation it is said that if boundary condition is not specified, it is automatically assumed to be Neumann 0). Therefore, I try to calculate the flow of the function at a random point of the boundary at the lower boundary:

D[sol[1, 1, z], z] /. z -> -10

And obtain -0.00125585, and not 0 as I expected. What is wrong with my solution?

---

Also, I am confused with the result not being smooth. See the following:

Plot[Evaluate[sol[1, 1, z], {z, -10, 0}]]

Another check shows that probably accuracy of the solution is unsufficient even at the boundary, however, I am not sure how to improve it without significant cost to the speed of computation:

Plot[{Evaluate[sol[1, y, 0]], y}, {y, 0, 100}, PlotRange -> All]

---

I tried using right-hand side in the following form: NeumannValue[- x y, z == 0] + NeumannValue[0, z != 0], which doesn't affect the answer.

Answer 1

I think this is somewhat an artefact of the discretization. Note that Neumann boundary conditions are only applied implicitely in the weak formulation of the PDE. In the weak formulation, the Neumann condition may be indeed fulfilled, but not in the way you compute it by differentiating with (interpolated) solution of the discrete system. Note, that we can decrease the mesh size for the discretization, e.g. as follows.

reg0 = Parallelepiped[{0, 0, -10}, {{100, 0, 0}, {0, 100, 0}, {0, 0, 10}}];
reg = DiscretizeRegion[reg0, MaxCellMeasure -> 5]; sol = 
 NDSolveValue[
  Laplacian[mu[x, y, z], {x, y, z}] - mu[x, y, z] == 
   NeumannValue[-x y, z == 0], mu, {x, y, z} \[Element] reg];

When I compute the Neumann differential as you did, then its absolute is smaller than before:

D[sol[1, 1, z], z] /. z -> -10

-0.000587538

This is somewhat an indicator for the fact that the two notions of Neumann boundary conditions (weakly and strongly formulated) coincide only in the limit $\text{mesh size} \to 0$.

Answer 2

I tried to solve the problem inside a smaller parallelepiped(based on the approach from Henrik Schumacher) because I assumed a scaling problem.

reg0 = Parallelepiped[{0,0, -1 }, {{1 , 0, 0}, {0, 1 , 0}, {0, 0, 1 }}];
reg = DiscretizeRegion[reg0 , MaxCellMeasure -> {"Volume" -> .001}] ;
sol = NDSolveValue[Laplacian[mu[x, y, z], {x, y, z}] - mu[x, y, z] == 
NeumannValue[x y, z == 0] + NeumannValue[0, z != 0],mu, {x, y, z} \[Element] reg]

The neumann conditions can be visualized

Plot3D[Evaluate[Derivative[0, 0, 1][sol][x, y, 0]], {x, 0, 1}, {y, 0,1}, PlotRange -> {-1, 1}] (* z=0: mu= x y *)

Plot3D[Evaluate[Derivative[0, 0, 1][sol][x, y, -1]], {x, 0, 1}, {y, 0,1}, PlotRange -> {-1, 1}]

and both plots agree well with the conditions given in NDSolve. Finally the examplary verification

Plot [ sol[1, 1, z] , {z, -1, 0 } , PlotRange -> {-1, 0}]

becomes smooth.

I'm still unsure what happens, perhaps scaling or meshing issue?

Answer 3

Up front, I can not explain exactly what is going on but I do not think this is a bug. Let's look at a simplified 2D version of this:

Needs["NDSolve`FEM`"]
w = 100;
l = 10;
r = Rectangle[{0, -l}, {w, 0}];
m = ToElementMesh[r,(*"MeshElementType"\[Rule]TriangleElement,*)
   "MeshOrder" -> 2, MaxCellMeasure -> 10];
(*m["Wireframe"]*)

nval = y;
sol = NDSolveValue[{Inactive[
        Div][-1. Inactive[Grad][mu[y, z], {y, z}], {y, z}] + 
      mu[y, z] == NeumannValue[nval, z == 0], 
    DirichletCondition[mu[y, z] == 0, z == -l]}, 
   mu, {y, z} \[Element] m];

I have made a few changes to the PDE. I put this in inactive form (does not make a difference) and I also added a DirichletCondition at z==0 - does not matter either but rules out that what we see is an issue due to an neumann only BC. It does not seem so.

D[sol[1, z], z] /. z -> -l
0.00008033698364150563`

Looks reasonable and not too far off compared to a Neumann zero value at that end.

Plot[sol[1, z], {z, -l, 0}, PlotRange -> All]

Note, that we also see this slight kink; but much less than with TriangleElement. This is also true for the 3D version. Using Cuboid and HexahedronElement makes this behave better.

Plot[Evaluate[{sol[y, z] - nval} /. {mu -> sol, z -> 0}], {y, 0, w}, 
 PlotRange -> All]

Note the drift in the difference between the neumann value and the computed value and also the jump at both ends. The jumps at both ends, I think, are there because the neumann zero BC also needs to be satisfied. The drift can be reduced is the mesh is refined but it remains.

Also, the

ContourPlot[sol[y, z], {y, 0, w}, {z, -l, 0}, PlotRange -> All]
Plot3D[sol[y, z], {y, 0, w}, {z, -l, 0}, PlotRange -> All]

look much smoother when using QuadElement (or HexahedronElement in 3D)

I have checked with other FEM software and I get a similar enough result to believe that this is not an implementation bug.