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]}
"SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}
$\endgroup$MeshCellMeasure
as a solution to your issue? $\endgroup$