Numerical solution of PDE does not satisfy initial conditions

Question

I want to solve the PDE for the variable $u(r,z,t)$

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

subject to the initial conditions $u(r,z,0) = re^{-(r^2+z^2)}$, $\frac{\partial u}{\partial t}(r,z,0) = 0$ and the boundary condition $u(0,z,t) = 0$.

The function $a(r)$ has the form of a step function.

This PDE describes an MHD wave propagating along an infinite cylinder embedded in an infinite medium, hence the physical domain is $0\leq r\leq\infty$, $-\infty\leq z\leq\infty$. Of course, I want to use a smaller numerical domain, e.g., $0\leq r\leq 4$, $-4\leq z\leq 4$. Although in this example I consider $0\leq t\leq 2$, where $t$ is time, I would like to solve the PDE for much longer times.

The code below is one of a few attempts I have tried to solve this problem with NDSolveValue and NDSolve. When I plot the numerically obtained $u$ for fixed $r$ and $z$ as a function of $t$, I obtain solutions that I know to be wrong, so my guess is that I am not correctly using NDSolveValue and NDSolve.

In the code below I simply plot the numerical solution and its error at $t=0$ along the line $r=r0$ and the line $z=z0$. There is a non-negligible difference between the numerical solution and the initial condition. In addition, I am quite surprised to find that this error changes if I modify the numerical domain in which the PDE is solved.

I will thank any ideas on how to improve the accuracy of the solution.

PDE = Derivative[0, 0, 2][u][r, z, t] == 
   a[r] (Derivative[2, 0, 0][u][r, z, t] + 
      1/r Derivative[1, 0, 0][u][r, z, t] - 1/r^2 u[r, z, t] + 
      Derivative[0, 2, 0][u][r, z, t]);
f[r_, z_] = r Exp[-(r^2 + z^2)];
ic = {u[r, z, 0] == f[r, z], Derivative[0, 0, 1][u][r, z, 0] == 0};
bc = u[0, z, t] == 0;

a[r_] = UnitStep[r - 1] + 1;

rmin = 0;
rmax = 4;
zmin = -4;
zmax = -zmin;
tmax = 2;

usol = NDSolveValue[{PDE, ic, bc}, 
   u, {r, rmin, rmax}, {z, zmin, zmax}, {t, 0, tmax},  
   Method -> {"MethodOfLines", "TemporalVariable" -> t, 
     "SpatialDiscretization" -> "FiniteElement"}, 
   InterpolationOrder -> All];

r0 = 1;
z0 = 0;
{Plot[{usol[r0, z, 0], f[r0, z]}, {z, -5, 5}, AxesLabel -> {"z", "u"}],
 Plot[f[r0, z] - usol[r0, z, 0], {z, -5, 5}, 
  AxesLabel -> {"z", "Error"}, PlotRange -> All]}
{Plot[{usol[r, z0, 0], f[r, z0]}, {r, 0, 4}, 
  AxesLabel -> {"r", "u"}],
 Plot[f[r, z0] - usol[r, z0, 0], {r, 0, 4}, 
  AxesLabel -> {"r", "Error"}, PlotRange -> All]}

Error at t=0 for r=0 using original code

Error at t=0 for r=0 using "MaxCellMeasure" -> 0.01

Answer 1

One way to improve this is to use a finer mesh via an option like:

"SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}

Since you have a discontinuity at r==1 you could also add

{"MaxCellMeasure" -> 0.01, "IncludePoints"->{{1,0}}}

To force a line in the mesh at r==1 Play a bit with that and see if it helps.