Neumann b.c are not satisfied in NDSolve PDE

Question

I am solving a PDE (p[x,t]) with 1 spatial dimension {x,-10,10} and 1 time dimension {t,0,100}. There is an external function, U. The code reads

U[x_] := Sin[ 12*2 Pi/20 x];

sol = NDSolve[{
     D[p[x, t], t] - D[p[x, t], x, x] + D[ U[x] p[x, t], x] - 1 + 
       p[x, t] == NeumannValue[0, x == -10] + NeumannValue[0, x == 10],
     p[x, 0] == 1
     }
    , p, {t, 0, 100}, {x, -10, 10}][[1]];

I want to impose dp/dx = 0 at x=-10 and x=10 at all times. Checking with:

state = NDSolve`ProcessEquations[{
     D[p[x, t], t] - D[p[x, t], x, x] + D[ U[x] p[x, t], x] - 1 + 
       p[x, t] == NeumannValue[0, x == -10] + NeumannValue[0, x == 10],
     p[x, 0] == 1
     }
    , p, {t, 0, 100}, {x, -10, 10}][[1]];

state["FiniteElementData"][
  "PDECoefficientData"]["DiffusionCoefficients"]
state["FiniteElementData"][
  "PDECoefficientData"]["ConservativeConvectionCoefficients"]
state["FiniteElementData"][
  "PDECoefficientData"]["LoadDerivativeCoefficients"]
(*
{{{{-1}}}}
{{{{0}}}}
{{{{0}}}}
*)

confirms that the Neumann value of 0 should only apply to the diffusive term, as the "ConservativeConvectionCoefficients" and "LoadDerivativeCoefficients" are 0. So the problem should be correctly imposing dp/dx = 0 at x=-10 and x=10. However

Table[D[(p /. sol)[x, t], x] /. {x -> -10}, {t, 0, 100}]
Table[D[(p /. sol)[x, t], x] /. {x -> 10}, {t, 0, 100}]

shows that the solution does not give a derivative of 0 for dp/dx at any of the timepoints for x=10,-10. Moreover the solution is just weirdly asymmetric, i.e. try:

Plot[(p /. sol)[x, 100], {x, -10, 10}]

How can I integrate this equation better and force the Neumann b.c to be satisfied at all the time points? Note: If I switch to:

U[x_] := Sin[2 Pi/20 x];

the dp/dx at the boundaries get closer to 0, like 10^-3, but this is not very good, and eventually, I want U to be more oscillatory than Sin[2 Pi/20 x].

Answer 1

You need to make the grid denser. Try e.g.

molfem[measure_ : Automatic] := {"MethodOfLines", 
      "SpatialDiscretization" -> {"FiniteElement", 
          "MeshOptions" -> MaxCellMeasure -> measure}};

sol = NDSolve[{D[p[x, t], t] - D[p[x, t], x, x] + D[U[x] p[x, t], x] - 1 + p[x, t] == 0, 
               p[x, 0] == 1}, 
    p, {t, 0, 100}, {x, -10, 10}, Method -> molfem[0.01]][[1]];

Table[D[(p /. sol)[x, t], x] /. {x -> -10}, {t, 0, 100}]
(* {0., 2.78726*10^-7, 2.78725*10^-7, 2.78725*10^-7, 2.78725*10^-7, ... *)

The NeumannValue[0, …] can be omitted, because they're the default setting.