# 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!

• For 2) Have a look at the message ref page: ref/message/NDSolve/femper. For 3) have a look at the FEMDocumentation Jun 1, 2020 at 9:54

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]}


• Thanks so much for the answer! It does indeed work. Just curious, if you have the time, why the other methods fail? Seems Inactive[Grad] not supporting not supporting the optional chart argument is just a bug but is there a way to solve a vector PDE in general or must it always be converted scalar form, like you did by taking Div of both sides? And further must I always convert it into an operator equation instead of a system of eqs in terms of the partials? DSolve handles the latter fine while NDSolve complains it's overdetermined. Thanks again! Jun 1, 2020 at 17:15
• @bRost03 Numerical methods are very sensitive, and as a rule it is better to turn system of equations to the classical form so that system can recognize it and apply specific algorithm (since we use NDSolve and not our own hand-made method). Jun 1, 2020 at 18:48