I'm trying to solve the set of equations below describing the flow of a pot of water being heated slightly. The equations are 2D axisymmetric in nature.
x0 = 1;
y0 = 1;
\[CapitalOmega] = Rectangle[{0, 0}, {x0, y0}];
\[Mu]water = 8.9*10^-4;
\[Alpha]water = 2.1* 10^-4;
kwater = 0.6;
Cwater = 4184;
\[Rho]water = 1000;
g = 9.78;
eqn1=D[vr[r,z,t],t]+vz[r,z,t]*D[vr[r,z,t],z]+vr[r,z,t]*D[vr[r,z,t],r]+D[p[r,z,t],r]/\[Rho]water-\[Mu]water (D[vr[r,z,t],z,z]+(D[r*D[vr[r,z,t],r],r]/r)-vr[r,z,t]/r^2)/\[Rho]water==0;
eqn2=D[vz[r,z,t],t]+vz[r,z,t]*D[vz[r,z,t],z]+vr[r,z,t]*D[vz[r,z,t],r]+D[p[r,z,t],z]/\[Rho]water-\[Mu]water (D[vz[r,z,t],z,z]+(D[r*D[vz[r,z,t],r],r]/r))/\[Rho]water==g \[Alpha]water T[x, y, t];
eqn3=D[vz[r,z,t],z]+D[r vr[r,z,t],r]/r==0;
eqn4=D[T[r,z,t],t]-kwater ((D[r*D[T[r,z,t],r],r]/r)+D[T[r,z,t],z,z])/(\[Rho]water Cwater)+vz[r,z,t] D[T[r,z,t],z]+vr[r,z,t] D[T[r,z,t],r]==0;
wall = DirichletCondition[{vr[r, z, t] == 0, vz[r, z, t] == 0}, r == 1];
reference = DirichletCondition[p[r, z, t] == 0, r == 0 && z == 0];
temperatures = {DirichletCondition[T[r, z, t] == 1, r == 1]};
bcs = {wall, reference, temperatures};
ic = {vz[r, z, 0] == 0, vr[r, z, 0] == 0, p[r, z, 0] == 0, T[r, z, 0] == 0};
Monitor[AbsoluteTiming[{xVel, yVel, pressure, temperature} = NDSolveValue[{eqn1, eqn2, eqn3,eqn4, bcs, ic}, {vz, vr, p, T}, {r, z} \[Element] \[CapitalOmega], {t, 0, 200},Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement","MeshOptions" -> {"MaxCellMeasure" -> 0.0000625}, "InterpolationOrder" -> {vz -> 2, vr -> 2, p -> 1, T -> 2}}}}, EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])];],currentTime]
After running the code, I get the error that Delay PDEs are not supported by mathematica. When I search online, most of the methods to fix this error is to use the method of lines which I have used in NDSolve, but still have not been able to fix this. Anyone knows what's wrong?