# Unable to solve Delay PDEs Error in Boussinesq Approximation

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?

• @user21 I think this is another parsing bug? Mar 17 at 2:51
• @xzczd I think the firs thing to do is to put it in inactive form. Mar 17 at 7:02
• @user21 is it possible to put it into inactive form? I don’t know if I can specify grad and div in terms of cylindrical coordinates while using the axisymmetric condition. Also out of curiosity what’s the difference between the inactive and active forms when it comes to solving? Mar 17 at 7:08
• Have a look at the DiffusionPDETerm (et al) ref pages, that have a "RegionSymmetry" -> "Axisymmetric" parameter. I have not tried this myself, you might need to fiddle a little with the equation too. Concerning the the inactive form, this specifies what the meaning of NeumannValue is, but also a few other things.... Mar 17 at 7:46
• ...See here and here, Not quite what you are looking for but perhaps also useful is the axisymmetric Tubular Reactor example. Mar 17 at 7:46