# How to extract particle velocity from acoustic pressure in an adjustable region?

I have set up an adjustable region in which to compute acoustic pressure, and would like to also know the particle velocity components, which are related to the gradient of the pressure. But I cannot figure out how to perform a gradient on the pressure function produced by ParametricNDSolveValue.

My region is axisymmetric, in cylindrical coordinates (r,z), between a top radiating surface and a lower reflecting paraboloid, with a reflecting spherical particle between. I seem to get plausible pressure using the following.

Needs["NDSolveFEM"];
rc = 0.03; (* Cylinder radius *)
(* hc = Cylinder height *)
(* fL = Parabolic reflector focal length *)
(* f = Frequency of radiating transducer *)
ca = 343.344; (*Sound speed in air *)
\[Rho]a = 1.20458 (*Air density *);
vars = {p[r, z], \[Omega], {r, z}};
pars = <|"MassDensity" -> \[Rho]a, "SoundSpeed" -> ca,
"RegionSymmetry" -> "Axisymmetric"|>;
air[hc_, fL_] =
ImplicitRegion[0 <= r <= rc && r^2/4/fL - fL <= z <= hc, {r, z}];
particle[rp_] =
ImplicitRegion[
r^2 + z^2 <= rp^2 && -rp <= z <= rp && r >= 0, {r, z}];
region[hc_, fL_, rp_] :=
DiscretizeRegion[RegionDifference[air[hc, fL], particle[rp]],
PrecisionGoal -> 7., AccuracyGoal -> Infinity];
mesh[hc_, fL_, rp_] := ToElementMesh[region[hc, fL, rp]];
pars, <|"SoundIncidentPressure" -> 1|>];
sideabsorb =
AcousticAbsorbingValue[
r == rc, {p[r, z], \[Omega], {r, z}}, <|
"AcousticSourceDistance" -> 1/(2*r), "SoundSpeed" -> ca,
"MassDensity" -> \[Rho]a|>];
pde[hc_] = {AcousticPDEComponent[vars, pars] ==
MaxGridSize[c_, \[Omega]_, Resolution_] :=
Module[{\[Lambda]}, \[Lambda] = 2 \[Pi] c/\[Omega];
\[Lambda]/Resolution];
Options[AcousticsOptions] = {"Resolution" -> 60};
AcousticsOptions[c_, \[Omega]_, opts : OptionsPattern[]] :=
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> {"Length" ->
MaxGridSize[c, \[Omega], OptionValue["Resolution"]]}}}};
pfun[hc_, fL_, rp_] :=
ParametricNDSolveValue[pde[hc],
p, {r, z} \[Element] region[hc, fL, rp], {\[Omega]},
AcousticsOptions[c, \[Omega]]];
Manipulate[
GraphicsGrid[{{
mesh[hc, fL, rp]["Wireframe"],
DensityPlot[
Evaluate[Re[pfun[hc, fL, rp][f*2*Pi][r, z]]], {r, 0,
rc}, {z, -fL, hc}, AspectRatio -> 2, PlotLegends -> Automatic],
DensityPlot[
Evaluate[Re[pfun[hc, fL, rp][f*2*Pi][r, z]]], {r, 0,
rc/2}, {z, -2*rp, 2*rp}, AspectRatio -> 2,
PlotLegends -> Automatic]
}}, ImageSize -> 600
], {{f, 22000, "frequency (Hz)"}, 20000,
30000}, {{hc, 0.06, "Height"}, 0.01,
0.1}, {{fL, 0.03, "Focal length"}, 0.01,
0.05}, {{rp, 0.003, "Particle radius"}, 0.001, 0.009}]

I would also like to see and work with the particle velocity components. But taking the Grad of the pressure doesn't give something I can plot - I get error statements to the effect that "True cannot be used as a variable", etc. So it seems that this is not the right way to get velocity:

Grad[Re[pfun[hc, fL, rp][f*2*Pi][r, z]], {r, z}];

I guess I don't understand how to work with the results of pfun. Experimenting, I can find this particular case without adjustable parameters where Grad works on pfun:

But what is the right way to extract velocity information within my adjustable region in the general case?

I did not find a way to use Grad on pfun, but I learned by experimentation that ParametricNDSolveValue, along with solving for the pressure p[r,z], could be asked to solve for its partial derivatives with respect to the coordinates, which would be the components of the gradient. For example,

dpdr = ParametricNDSolveValue[pde,
D[p[r, z], r], {r, z} \[Element] region, {\[Omega]},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> {"Length" ->
maxgridsize}}}}];

I'm not sure where this capability of ParametricNDSolveValue is mentioned in documentation, but I did not find any examples, so it surprised and delighted me.

From these partial derivatives I can get the components of particle velocity in the radial and axial directions, vr and vz. In a simplified Manipulate example below they are combined to plot the magnitude of the velocity:

Needs["NDSolveFEM"];
rc = 0.03; (* Cylinder radius *)
(* hc = Cylinder height *)
(* fL = Parabolic reflector focal length *)
(* f = Frequency of radiating transducer *)
ca = 343.344; (*Sound speed in air *)
\[Rho]a = 1.20458 (*Air density *);
vars = {p[r, z], \[Omega], {r, z}};
pars = <|"MassDensity" -> \[Rho]a, "SoundSpeed" -> ca,
"RegionSymmetry" -> "Axisymmetric"|>;
Manipulate[
{air :=
ImplicitRegion[0 <= r <= rc && r^2/4/fL - fL <= z <= hc, {r, z}];
particle :=
ImplicitRegion[
r^2 + z^2 <= rp^2 && -rp <= z <= rp && r >= 0, {r, z}];
region :=
DiscretizeRegion[RegionDifference[air, particle],
PrecisionGoal -> 5., AccuracyGoal -> Infinity];
pars, <|"SoundIncidentPressure" -> 1|>];
sideabsorb :=
AcousticAbsorbingValue[
r == rc, {p[r, z], \[Omega], {r, z}}, <|
"AcousticSourceDistance" -> 1/(2*r), "SoundSpeed" -> ca,
"MassDensity" -> \[Rho]a|>];
pde := {AcousticPDEComponent[vars, pars] == topradiate + sideabsorb};
maxgridsize := ca/f/12;
\[Omega]0 := f*2*Pi; (*Angular frequency*)
pfun =
ParametricNDSolveValue[pde,
p, {r, z} \[Element] region, {\[Omega]},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> {"Length" ->
maxgridsize}}}}];
dpdr =
ParametricNDSolveValue[pde,
D[p[r, z], r], {r, z} \[Element] region, {\[Omega]},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> {"Length" ->
maxgridsize}}}}];
dpdz =
ParametricNDSolveValue[pde,
D[p[r, z], z], {r, z} \[Element] region, {\[Omega]},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> {"Length" ->
maxgridsize}}}}];
vr = -I/(\[Rho]a*\[Omega]0)*dpdr[\[Omega]0];
vz = -I/(\[Rho]a*\[Omega]0)*dpdz[\[Omega]0];
GraphicsRow[{
DensityPlot[
Evaluate[Norm[pfun[\[Omega]0][r, z]]], {r, 0, rc}, {z, -fL,
hc}],
DensityPlot[
Evaluate[Norm[{vr, vz}]], {r, 0, rp*3/2}, {z, -rp*3/2, rp*3/2}]
}]},
{{f, 22000, "Frequency"}, 20000, 30000,
Appearance -> {"Labeled"}}, {{hc, 0.06, "Height"}, 0.01, 0.1,
Appearance -> {"Labeled"}}, {{fL, 0.03, "Focal length"}, 0.01, 0.06,
Appearance -> {"Labeled"}}, {{rp, 0.005, "Particle radius"}, 0.001,
0.009, Appearance -> {"Labeled"}}]