Solve 2D diffusion equation in polar coordinates

Question

I am trying to solve the diffusion equation in polar coordinates:

$\frac{\partial u}{\partial t} = D \left( \frac{\partial^2 u}{\partial r^2} + \frac{1}{r} \frac{\partial u}{\partial r}\right)$

subject to there being a Dirichlet boundary condition at $r=10$, and there being no-flux (Neumann) at $r=0$. I am starting with $u(0,r)=e^{-10r}$, although would eventually like to move to a delta function at the origin.

I have tried the following (where I have set $D=0.01$):

d = 0.01;

op2Diff = -d*(Derivative[0, 1][u][t, r]/r + Derivative[0, 2][u][t, r]) + 
    Derivative[1, 0][u][t, r];

uSoln2Diff = 
  NDSolveValue[{op2Diff == 0 + NeumannValue[0, r == 0], 
DirichletCondition[u[t, r] == 0, r == 10], u[0, r] == Exp[-10 r]},
u, {t, 0, 1000}, {r, 0, 10}];

I get the following error:

    NDSolveValue::femcnmd: "The PDE coefficient {{1.,-(0.01/r)}} does not evaluate to a numeric matrix of dimensions {1,2}."

Does anyone know what this means, and how to rectify the issue?

Also, following my point above - does anyone know how to specify a delta function at the origin as an initial condition? I have tried to specify $u(0,r)=DiracDelta[r]$, but this results in an error saying that there isn't a numeric derivative at this time point (even for simpler PDE examples than the one here, which work with other initial conditions).

Sorry for the poor formatting of the Mathematica code above - I don't know how to make it look neat (please let me know if you do).

Best,

Ben

Answer 1

The error is occurring, because Derivative[0, 1][u][t, r]/r evaluates to 0/0 at r = 0. A simple way to handle this is to move the inner boundary to some small value, say, r0 = 1/1000. With this change and a few simplifications,

uSoln2Diff = NDSolveValue[{op2Diff == 0, (D[u[t, r], r] /. r -> r0) == 0,
    u[t, 10] == 0, u[0, r] == Exp[-10 r] - Exp[-100]}, u, {t, 0, 1000}, {r, r0, 10}];
Plot3D[uSoln2Diff[t, r], {t, 0, 1000}, {r, r0, 10}, PlotRange -> All, 
    AxesLabel -> {t, r, u}]

Note that u[0, r] has been modified infinitesimally, and perhaps unnecessarily, so that it is consistent with u[t, 10] == 0. There still may be an inconsistency between initial and boundary conditions at r = r0, but the answer looks reasonable.

Answer 2

This works out of the box in version 11.3 (possibly also earlier versions):

d = 0.01;

op2Diff = -d*(Derivative[0, 1][u][t, r]/r + Derivative[0, 2][u][t, r]) + 
    Derivative[1, 0][u][t, r];

uSoln2Diff = 
  NDSolveValue[{op2Diff == 0 + NeumannValue[0, r == 0], 
DirichletCondition[u[t, r] == 0, r == 10], u[0, r] == Exp[-10 r]},
u, {t, 0, 1000}, {r, 0, 10}];

The message you saw was caused by the verification of the PDE coefficients that where parsed. This verification happened at the coordinate 0, which in this case caused the message and the rejection of the coefficient. In newer versions this verification happens with a coordinate from some place inside the domain. The actual integration did not need to evaluate the coefficient at the singularity.

It is possible to change the coordinate that is used for the verification of the coefficient. This can be done with specifying VerificationData:

NDSolveValue[{op2Diff == 0 + NeumannValue[0, r == 0], 
  DirichletCondition[u[t, r] == 0, r == 10], 
  u[0, r] == Exp[-10 r]}, u, {t, 0, 1000}, {r, 0, 10}, 
 Method -> {"MethodOfLines", 
   "SpatialDiscretization" -> {"FiniteElement", 
     "InitializePDECoefficientsOptions" -> {"VerificationData" -> \
{"Coordinate" -> {0.}}}}}]

Will give the message again. VerificationData is documented in the options section of InitializePDECoefficients.