Inconsistent boundary and initial conditions: BC ignored altogether

Question

Consider the following diffusion-decay equation with von Neumann b/c in the origin and Dirichlet at the other boundary:

pde = D[u[x, t], t] == d*D[u[x, t], x, x] - \[Delta]*u[x, t]
bc = {D[u[x, t], x] == -j0 /. x -> 0, u[x, t] == 0 /. x -> length}
ic = u[x, t] == u0 /. t -> 0
params = {d -> 10, \[Delta] -> 1, length -> 600, u0 -> 200, j0 -> 4000}
tFinal = 150;
soln = NDSolve[Flatten@{pde, bc, ic} /. params, u, {t, 0, tFinal}, {x, 0, length /. params}]

Running NDSolve as above generates this Warning:

NDSolve::ibcinc: Warning: Boundary and initial conditions are inconsistent. 

I understand that the boundary and initial conditions are inconsistent. To solve this problem I looked through this webpage:

https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html#1306392612

I noticed, however, that changing j0 had no effect on the attained numerical solution and, in fact, leaving the symbol j0 undefined does not lead NDSolve to complain about non-numerical derivatives and such. This leads me to believe that NDSolve ignores somehow the boundary condition in the origin altogether, making it somewhat hard to debug the warning message about inconsistent boundary and initial conditions.

To clarify my point about j0: I don't understand how NDSolve manages to ignore the boundary condition in the origin altogether without running into trouble. If I comment out the left b/c, NDSolve complains:

NDSolve::bcart: Warning: An insufficient number of boundary conditions have been specified for the direction of independent variable x. Artificial boundary effects may be present in the solution.

Edit: One solution I found is to turn j0 into a Michaelis-Menten type of function of time. Here, I also incorporate a valid point made in the comments: the Dirichlet b/c on the right of the domain was also inconsistent with the initial condition.

pde = D[u[x, t], t] == d*D[u[x, t], x, x] - \[Delta]*u[x, t]
bc = {D[u[x, t], x] == -((j0*t)/(1 + t)) /. x -> 0, u[x, t] == 0 /. x -> length}
ic = u[x, t] == u0 /. t -> 0
params = {d -> 100, \[Delta] -> 1, length -> 600, u0 -> 0, j0 -> 40}
tFinal = 150;
soln = NDSolve[Flatten@{pde, bc, ic} /. params, u, {t, 0, tFinal}, {x, 0, length /. params}, Method -> {"MethodOfLines", 
"SpatialDiscretization" -> {"TensorProductGrid", 
  "MinPoints" -> 1000}}]

Your help and suggestions are greatly appreciated!

Answer 1

The solution has actually been implied on the page you linked. (Well, I admit that tutorial isn't that readable, I myself didn't notice the solution therein until recently.) To fix the problem, you need to add option Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}} to NDSolve, because NDSolve will automatically set "ScaleFactor" -> 0 for non-Dirichlet boundary condition, when i.c. and b.c. are inconsistent, this will make the non-Dirichlet inconsistent b.c. be severely ignored. (In your case since the right hand side of your Neumann condition D[u[x, t], x] == -j0 /. x -> 0 isn't a function of t, it's completely ignored. ) For more information you can revisit the tutorial.

After adding the option, the solution is obviously influenced by j0:

pde = D[u[x, t], t] == d*D[u[x, t], x, x] - δ*u[x, t];
bc = {D[u[x, t], x] == -j0 /. x -> 0, u[x, t] == 0 /. x -> length};
ic = u[x, t] == u0 /. t -> 0;
params = {d -> 10, δ -> 1, length -> 600, u0 -> 200};
tFinal = 150;
solp = ParametricNDSolveValue[Flatten@{pde, bc, ic} /. params, 
  u, {t, 0, tFinal}, {x, 0, length /. params}, j0, 
  Method -> {"MethodOfLines", 
    "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}]
Plot[solp[#][x, 100] & /@ Range[4] // Evaluate, {x, 0, length/10 /. params}, 
 PlotRange -> All]