NDSolve PDE trouble with overdetermined, inactive, inconsistent equation dimensions

Question

TL;DR

How can I solve $\left(\vec{v}\cdot\nabla\right)\vec{v}=-g\hat{z}-\nabla p$ numerically for $p(r,\theta,z)$ given an explicit $\vec{v}(r,z)$ and $p(r_0,\theta,z)=g z$?

Issue

I have a PDE which I would like to solve but I am having a lot of difficulty. It's a relatively simple looking equation I'd like to solve for $p$ $$ \left(\vec{v}\cdot\nabla\right)\vec{v}=-g\hat{z}-\nabla p $$ $g$ is a constant while $\vec{v}$ and $p$ are best handled in cylindrical coordinates, $\{r,\theta,z\}$, as $\partial_\theta p=\partial_\theta \vec{v}=0$ and $\vec{v}\cdot\hat{r}=\vec{v}\cdot\hat{z}=0$. I know that for a given $r_0$ I have $p(r_0,\theta,z)=g z$. I have two different $\vec{v}$s I'd like to solve this for. One is quite simple and can be done analytically with DSolve, this worked well. The other $\vec{v}$ is an infinite sum of Bessel functions which Mathematica cannot handle so I want to solve this numerically.

I have attempted to do this in operator form as well as in component form but I cannot seem to get Mathematica to do it. Here's an example $\vec{v}$, the actual expression is much uglier.

v[r_, z_, nmax_] := {0, Sum[(BesselI[1, (2 n - 1) r] + BesselK[1, (2 n - 1) r]) Sin[(2 n - 1) z], {n, 1, nmax}], 0};

Operator method:

If I try to solve this in operator form I have problems.

diffEq[nmax_] := N@Block[{g = {0,0,1}, 
diff = (v[r, z, nmax].Inactive[Grad][#, {r, t, z},"Cylindrical"]) &@v[r, z, nmax]}, 
   Inactive[Grad][p[r, t, z], {r, t, z},"Cylindrical"] == (-diff - g)]
NDSolve[{diffEq[2], DirichletCondition[p[r, t, z] == z, r == 2]}, p, {r, 1, 2}, {z, 1, 2}, {t, 1, 2}]
bc[r0_] = DirichletCondition[p[r, t, z] == z, r == r0];

Gives me an error

Inactive::argrx: Inactive[Grad] called with 3 arguments; 2 arguments are expected.

Which shouldn't be right, this seems like a problem with Mathematica. But no matter, I don't really care about $\theta$ and the Cartesian gradient works the same for $r$ and $z$ so we try

diffEq[nmax_] := N@Block[{g = {0,0,1}, 
diff = (v[r, z, nmax].Inactive[Grad][#, {r, t, z}]) &@v[r, z, nmax]}, 
   Inactive[Grad][p[r, t, z], {r, t, z}] == (-diff - g)]
NDSolve[{diffEq[2], bc[2]}, p, {r, 1, 2}, {z, 1, 2}, {t, 1, 2}]

And I get

NDSolve::femper: PDE parsing error of ... Inconsistent equation dimensions.

Which seems to be Mathematica complaining that the LHS is a single operator expression while the RHS is a list because $g\hat{z}$ is a vector? This happens whether or not I include the boundary condition. I found this question which wasn't very helpful. How can this be overcome?

Component method:

I also tried to evaluate the gradients and solve the resulting system.

NDSolve[Flatten[{Activate@diffEq[1],bc[2]}], p, {r, 1, 2}, {z, 1, 2}, {t, 1, 2}]

gives

NDSolve::overdet: There are fewer dependent variables, {p[r,t,z]}, than equations, so the system is overdetermined.

Which kind of makes sense because I have three equations for 1 unknown, but it's really not overdetermined. This happens whether or not I include the boundary condition. I can drop the equation for $\theta$ (aka t) since it's trivial.

diffEq[nmax_] := Block[{g = {0, 0, 1},
  diff = (v[r, z, nmax].Grad[#, {r, t, z}]) &@v[r, z, nmax]}, 
  Grad[p[r, z], {r, t, z}] == (-diff - g)]
NDSolve[(Activate@diffEq[1])[[All, {1, 3}]], p, {r, 1, 2}, {z, 1, 2}]

Still throws the same error as expected but

DSolve[(Activate@diffEq[1])[[All, {1, 3}]], p, {r, 1, 2}, {z, 1, 2}]

Works fine and gives the solution (although this doesn't actually work because I need diffEq[n] where n is bigger than 1). So what gives? Why does NDSolve refuse to even try something that DSolve can handle without complaint?

Seeking:

Any or all of the following

1. Working code to solve my PDE
2. Insight into what I'm doing wrong
3. Tips, tricks or hacks when dealing with PDE's in Mathematica

Thanks!

Answer 1

There is no dependence on t in this problem. Also it could be easily turn to the second order equation. For first test the code should be like this

nmax = 2; v = {0, 
  Sum[(BesselI[1, (2 n - 1) r] + 
      BesselK[1, (2 n - 1) r]) Sin[(2 n - 1) z], {n, 1, nmax}], 0};

g = {0, 0, 1}; diff = 
 Div[(v.Grad[#, {r, \[Theta], z}, "Cylindrical"]) &@v, {r, \[Theta], 
   z}, "Cylindrical"]; eq = 
 Laplacian[p[r, z], {r, \[Theta], z}, "Cylindrical"] + diff; reg = 
 ImplicitRegion[1 <= r <= 2 && 1 <= z <= 2, {r, z}];
sol = NDSolveValue[{eq == 0, 
    DirichletCondition[p[r, z] == z, r == 2]}, 
   p, {r, z} \[Element] reg];

Visualization with boundary condition

{DensityPlot[sol[r, z], {r, z} \[Element] reg, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  PlotPoints -> 50, FrameLabel -> Automatic], 
 Plot[sol[2, z], {z, 1, 2}, AxesLabel -> Automatic]}