This is continuation from my previous post How to ensure for a solution of NDSolve to be positive? [https://mathematica.stackexchange.com/questions/278777/how-to-ensure-for-a-solution-of-ndsolve-to-be-positive]
I'am trying to solve Fokker-Plank kinetic equation for the test particles injected into a Maxwellian plasma. It is well known in plasma physics and has many applications. In particular, this equation describes the evolution of a beam of fast neutral atoms that are injected into a dense target plasma to heat it. Neutral atoms are ionized in the plasma and, after being converted into ions, are confined in the plasma by magnetic fields. My question will be about how to correctly formulate the boundary conditions for solving this equation using NDSolveValue. But in order to substantiate this question, it seems to me necessary to explain the physical formulation of the problem, at least in a minimal amount. Those who are willing to believe that the problem is set correctly can skip the next section.
Problem
Usually, the density of fast ions is small compared to the density of the target plasma, so that the scattering of fast ions by fast ions can be neglected, which exactly corresponds to the test particle approximation. The evolution of fast ion distribution is described by PDE: $$ \frac{\partial }{\partial t}f = -\frac{1}{v^{2}} \left[ \frac{\partial}{\partial v} \left( v^{2}h'(v) \left( f(t, v,\theta ) + (2\mu/9\pi)^{1/3}v_{c}^{2}\frac{1}{v}\frac{\partial}{\partial v}f(t, v,\theta) \right) \right) - \frac{v_c^3}{2v}g'(v) \frac{1}{\sin(\theta )}\frac{\partial}{\partial\theta}\left( \sin(\theta )\frac{\partial}{\partial\theta}f(t, v,\theta ) \right) \right] + S(v,\theta ) - R(v,\theta )f(t,v,\theta ) , $$ where $\mu=m_{e}/m_{i}$ (ratio of electron mass to ion mass), $v_{c}=(3\sqrt{\pi}\mu/4)^{1/3}v_ {e}$, $v_{e}$ is the thermal velocity of the target plasma electrons, $S(v,\theta )$ is the source of fast ions, $R(v,\theta )$ is the ion loss rate (only in the loss cone, that is, at $\sin(\theta )^{2}<1/M)$, where $M>1$ is the so-called mirror ratio.
In the simplest approximation, it is assumed that the velocity of fast ions falls within the interval $v_{i} \ll v \ll v_{e}$. Then the Rosenbluth potentials $h(v)$ and $g(v)$ are approximated by very simple formulas: $v^2h'(v)=v^{3}+v_{c}^{3}$, $g'(v)=1$. Under these conditions, the second term in parentheses is often discarded, thereby ignoring energy diffusion. It is also often assumed that fast ions, which, due to the angle scattering (described by the last term in square brackets), penetrate into the loss cone, are almost instantly lost (escape through the magnetic mirrors). Then the equation is solved without the last term, but with a zero-boundary condition on the margin between the loss cone and region of confinement at $\theta_{LC}=\arcsin(1/\sqrt{M})$: $ f(v,\theta_{LC})=0 $. Then we come to the equation $$ \frac{\partial }{\partial t}f = -\frac{1}{v^{2}} \left[ \frac{\partial}{\partial v} \left( (v^{3}+v_{c}^{3}) \left( f(t, v,\theta ) \right) \right) - \frac{v_{c}^{3}}{2v} \frac{1}{\sin(\theta )}\frac{\partial}{\partial\theta}\left( \sin(\theta )\frac{\partial}{\partial\theta}f(t, v,\theta ) \right) \right] + S(v,\theta ) . $$ Ultimately, we will be interested in the steady-state solution of the above equation, when the distribution function does not depend on time. We need the time derivative on the left side of the equation only to clarify the meaning of the terms on the right side. Integrating both parts of the equation over the volume of the spherical layer $v_{min}<v<v_{max}$, $-\pi + \theta _{LC} <\theta <\theta_{LC}$, on the left side we get the rate of changes in the total number of fast ions in this layer, $dN/dt$. On the right side, the integral of the divergence in space velocity is converted to an integral over the surface of the particle flux through the corresponding parts of the boundary of the spherical segment. So we can conclude that $$ \mbox{Flux}_{v} =(v^{3}+v_{c}^{3})f(t, v,\theta ) $$ has the meaning of a particle flow along the radius (that is, along $v$), and $$ \mbox{Flux}_{\theta } = - \frac{v_c^{3}}{2v} \left( \sin(\theta )\frac{\partial}{\partial\theta}f(t, v,\theta ) \right) $$ has the meaning of the flow of particles through the boundary of the loss cone.
If we neglect the angle scattering for a while, then it is easy to see that in the region of velocities $v<v_{0}$ is less than the velocity $v_{0}$ of injected atoms, the solution of the stationary kinetic equation is proportional to $$ f\propto \frac{1}{v^{3}+v_{c}^{4}} .$$ In other words, in this case we have a solution with a constant flow from the source of particles into the region $v=0$. But there is no ion absorber here, so this solution is meaningless. Accounting for angular scattering reduces the flux to zero in velocity $v$, but not completely.
It is important that the first term taken into account (with the derivative with respect to velocity) describes the deceleration of ions, therefore, by discarding the second, diffusion, term in the original equation, it is not possible to construct a solution with a flow into the region of infinite velocities. Thus, we conclude that discarding the second term, which in the original version of the equation describes energy diffusion, leads to an insoluble contradiction, at least in the region of low velocities, $v\leq v_{i}$.
To resolve this contradiction, one can try to (1) restore the diffusion term and/or (2) use more precise formulas for the Rosenbluth potentials, which were also true for $v<v_{i}$.
Our numerical experiments have previously shown that method (2) only exacerbates this contradiction. Therefore, in the next section, we describe our attempts to follow method (1) with a slight refinement of the formulas for the Rosenbluth potentials.
How to nullify a stream at zero?
The original equation has singularities at $v=0$ and $\sin(\theta)=0$. The second of them can be destroyed by passing to the variable $\xi=\cos(\theta )$. After multiplying all terms of the equation by $v^{2}$, we get $$ 0 = - \left[ \frac{\partial}{\partial v} \left( (v^{3}+v_{c}^{3}) \left( f(v,\xi) + (2\mu/9\pi)^{1/3}v_{c}^{2}\frac{1}{v}\frac{\partial}{\partial v}f(v,\xi) \right) - \left( \frac{v_c^3}{2 v} \right) \frac{\partial}{\partial\xi} \left( (1-\xi^{2}) \frac{\partial}{\partial\xi}f(v,\xi) \right) \right) \right] + v^{2}S(v,\xi) - v^{2}R(v,\xi)f(v,\xi) . $$
The NDSolve documentation recommends writing the equation in the inactive form in case of difficulty:
$$
\nabla.\left(
-c.\nabla F - \alpha F
\right)
+ a F
=
f
,$$
where $c$ are called DiffusionCoefficient, $\alpha $ are ConservativeConvectionCoefficients, $a$ ReactionCoefficients, and $f$ are LoadCoefficients. In the code below Diffusion Coefficient are given by
$$
Dei = \left(
\begin{array}{cc}
-\frac{\sqrt[3]{\frac{2}{\pi }} \sqrt[3]{\mu }
\left(u^3+1\right)}{3^{2/3} u} & 2
\text{D12}(u,\xi ) \\
-2 \text{D12}(u,\xi ) & \frac{\xi ^2-1}{2 u} \\
\end{array}
\right)
,
$$
Conservative Convection Coefficients are given by
$$
Aei
=
\left\{2 \text{D12}^{(0,1)}(u,\xi )-u^3-1,-2
\text{D12}^{(1,0)}(u,\xi )\right\}
,
$$
Reaction Coefficients come from rate of loss
$$
R = -\frac{\pi \sqrt{\frac{M-1}{M}} u^3 \theta
\left(\xi ^2+\frac{1}{M}-1\right)}{2 {\tau c} \log \left(\frac{16
(M-1)}{M \left(M \left(\xi
^2-1\right)+1\right)}\right)}
,$$
and Load Coefficients come from source of particles,
$$
S
=
\frac{S0 }{2 \pi\, \Delta\xi\, \Delta u}
e^{-\frac{(\xi -\xi0)^2}{2 \Delta \xi
^2}-\frac{(u-u0)^2}{2 \Delta u}^2}
.
$$
DependentVariables is the distribution function
$
F[u,\xi],
$$
where $u=v/v_{c}.$
A couple more explanations. When reconstructing the diffusion coefficients and conservative convection coefficients from the original equation, it turned out that they are determined up to an arbitrary function $D12[u,\xi]$. It is not included in the equation, but is included in the expression for the flow. In the code below, this function is set to zero. The second clarification refers to the $iExp[x]$ function, which was introduced to block an error that occurs when $Exp[x]$ function has large negative arguments.
We solve the kinetic equation twice. First, with an empty loss cone, which formally corresponds to an infinitely small time of particle escape from the cone, $\tau c=0$. However, division by zero in Reaction Coefficients $R$ would surprise the Mathematician, so we solve the equation only in the plasma confinement region $umin<u<umax$, $0<\xi<\sqrt{1-1/M}$, where $R=0$. On the boundaries of this region, we impose the Dirichlet boundary condition $F=0$ for $u=umax$ and for $\xi=\sqrt{1-1/M}$. The boundary $\xi=0$ corresponds to the plane of symmetry of the containment area, so here we set the Neumann boundary condition to zero, which should theoretically correspond to the vanishing of the angular projection of the flux vector
$$
Fei
= -Dei . \nabla_{u,\xi}F-Aei*F
.$$
At the boundary $u=umin$ near the zero velocity modulus, we also set the zero Neumann boundary condition, since we want to say that there is no sink of particles into the region of zero velocity. But this is precisely what is not obtained in our calculations. We denote the solution of this problem as $fvLC$. To build graphs, we also introduce the function $fsLC$, which complements $fvLC$ with zero in the area outside its definition.
For the second time, we solve the equation with a small but finite particle escape time $\tau c$. In this case, we are looking for a solution for all admissible values of the angular variable $\xi$ from $0$ to $1$, and at the boundary $\xi=1$ we set the zero Neumann boundary condition, since this is also a plane of symmetry in the space of velocities.
The calculation results are presented in the form of 3D graphs of the function itself, 2D graphs of the density of the function. These graphs look quite plausible at first glance. However, the magnitude of the peak near zero velocity is confusing. It is especially large in the graph on the right, where the result of the second calculation is presented.
The third pair of plots shows a 2D plot of the radial component of the flux vector. Finally, the fourth pair of graphs depicts this radial component as a function of the angular coordinate $\xi$ for several values of the velocity $u$ near $umin$. From this graph it is best seen that there is a huge particle flux at zero velocity, which I think, wrong.
The right graph in the fourth row convincingly demonstrates that, near zero velocity, the found ion distribution function is isotropic (does not depend on the angular variable), as it should be, but it is not equal to the Maxwellian one, which was the case when the radial flow was equal to zero.
What is wrong in my reasoning or in my code? Why does the radial flow in the direction of velocity zero not vanish at the inner boundary $u=umin$, but rather increase as it approaches this boundary? Maybe the boundary conditions here are formulated incorrectly?
Updated 2023-02-15
I plotted the radial flow versus speed at a fixed angle equal to the angle of the source. It can be seen from the graph that the flux decreases linearly from the maximum at the source to zero at zero velocity. However, near zero velocity, a numerical instability appears in the form of an oscillator. The amplitude and frequency of oscillations decreased after changing the variable u -> w=u^2. Thus, the question of how to correctly set the BC must be recognized as obsolete. The boundary condition works as expected. But the question arises how to prevent instability near u=0.
MWE
First, some preliminary definitions:
ClearAll[iExp]
iExp[x_] /; (x >= -676.) = Exp[x]
iExp[x_] /; (x < -676.) = 0
Derivative[1][iExp][x_] = iExp[x]
Derivative[2][iExp][x_] = iExp[x]
Second, define coefficients of equation:
D12Rule = {D12 -> Function[{u, \[Xi]}, 0]};
resetRule = {Exp[x_] -> iExp[x]};
Dei =
{{-(((2/\[Pi])^(1/3) (u^3 + \[Eta]) \[Mu]^(1/3))/(
3^(2/3) u)), (1 + \[Eta]) D12[
u, \[Xi]]}, {-((1 + \[Eta]) D12[
u, \[Xi]]), (\[Eta] (-1 + \[Xi]^2))/(2 u)}} /. D12Rule;
Aei =
{-u^3 - \[Eta] + (1 + \[Eta])
\!\(\*SuperscriptBox[\(D12\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[u, \[Xi]], -((1 + \[Eta])
\!\(\*SuperscriptBox[\(D12\),
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[u, \[Xi]])} /. D12Rule;
Fei[u_, \[Xi]_, {\[Mu]_, \[Eta]_}] = -Dei .
Grad[F[u, \[Xi]], {u, \[Xi]}] - Aei *F[u, \[Xi]];
gradF[u_, \[Xi]_, {\[Mu]_, \[Eta]_}] = -Grad[F[u, \[Xi]], {u, \[Xi]}];
inactiveFei[
u_, \[Xi]_, {\[Mu]_, \[Eta]_}] = -Dei .
Inactive[Grad][F[u, \[Xi]], {u, \[Xi]}] -
Inactive[Times][Aei, F[u, \[Xi]]] /. D12Rule;
inactiveEqn2ei$[
u_, \[Xi]_, {\[Mu]_, \[Eta]_, S0_, \[Tau]c_, M_,
u0_, \[CapitalDelta]u_, \[Xi]0_, \[CapitalDelta]\[Xi]_}] =
Inactive[
Div][-Dei . Inactive[Grad][F[u, \[Xi]], {u, \[Xi]}] -
Inactive[Times][Aei, F[u, \[Xi]]] /. D12Rule, {u, \[Xi]}] + (
E^(-((u - u0)^2/(2 \[CapitalDelta]u^2)) - (\[Xi] - \[Xi]0)^2/(
2 \[CapitalDelta]\[Xi]^2)) S0)/(
2 \[Pi] \[CapitalDelta]u \[CapitalDelta]\[Xi]) - (
Sqrt[(-1 + M)/
M] \[Pi] u^3 F[u, \[Xi]] UnitStep[-1 + 1/M + \[Xi]^2])/(
2 \[Tau]c Log[(16 (-1 + M))/(M (1 + M (-1 + \[Xi]^2)))]) /.
resetRule
activeEqn2ei$[
u_, \[Xi]_, {\[Mu]_, \[Eta]_, S0_, \[Tau]c_, M_,
u0_, \[CapitalDelta]u_, \[Xi]0_, \[CapitalDelta]\[Xi]_}] =
Div[-Dei . Grad[F[u, \[Xi]], {u, \[Xi]}] -
Times[Aei, F[u, \[Xi]]] /. D12Rule, {u, \[Xi]}] + (
E^(-((u - u0)^2/(2 \[CapitalDelta]u^2)) - (\[Xi] - \[Xi]0)^2/(
2 \[CapitalDelta]\[Xi]^2)) S0)/(
2 \[Pi] \[CapitalDelta]u \[CapitalDelta]\[Xi]) - (
Sqrt[(-1 + M)/
M] \[Pi] u^3 F[u, \[Xi]] UnitStep[-1 + 1/M + \[Xi]^2])/(
2 \[Tau]c Log[(16 (-1 + M))/(M (1 + M (-1 + \[Xi]^2)))]) /.
resetRule
Finally, solve equations and draw results:
ClearAll[F, fv, fs, fvLC, fsLC, fp, fpLC, \[Tau]c, \[Tau]s, M, \
\[Xi]0, u0, \[Xi]0]
With[{u0 = 5, umin = 0.002,
umax = 7, \[CapitalDelta]\[Xi] = 0.04, \[CapitalDelta]u = 0.08,
M = 8, \[Tau]c = 1./100, \[Tau]s = 1,
S0 = 1, \[Xi]0 = Cos[45. Degree]
, \[Mu] = 10*Quantity[1, "ElectronMass"]/Quantity[1, "ProtonMass"],
vc = 1, \[Eta] = 1
, eqn = inactiveEqn2ei$, flux = Fei(*, flux=gradF*)
, pltRange = Full
, pltRange3D = {-0.001, 1}, pltRangeDens = {-0.001, 1},
pltRange1D = Automatic, pltRange3DFlux =(*{-10,10}*)Automatic
, iSize = 280}
,
(*===========================================================================*)
(* Helper functions ===========================================================*)
getStateObject[eq_] :=
NDSolve`ProcessEquations[
eq, {u, \[Xi]}, {u, umin, umax}, {\[Xi], 0, Sqrt[1 - 1/M]}
(*,Method\[Rule]{"MethodOfLines",
"SpatialDiscretization"\[Rule]{"FiniteElement"}}*)
][[1]];
getEquations[eq_] := Module[{temp},
temp = getStateObject[eq];
temp = temp["FiniteElementData"];
temp = temp["PDECoefficientData"];
(# -> temp[#]) & /@ {"LoadCoefficients",
"LoadDerivativeCoefficients", "DiffusionCoefficients",
"ConservativeConvectionCoefficients", "ConvectionCoefficients",
"ReactionCoefficients", "DampingCoefficients"(*,
"MassCoefficients","Stationary"*)}];
(*============================================================================*)
(* Empty Loss Cone *)
eqn40LC = {(eqn[
u, \[Xi], {\[Mu], \[Eta], S0, \[Infinity], M,
u0, \[CapitalDelta]u, \[Xi]0, \[CapitalDelta]\[Xi]}] /. {Exp[
x_] -> iExp[x]}) ==
0 + NeumannValue[0, (\[Xi] == 0) && (umin <= u <= umax)] +
NeumannValue[0, (0 <= \[Xi] <= Sqrt[1 - 1/M]) && (u == umin)]
, DirichletCondition[F[u, \[Xi]] == 0.,
u == umax && 0 <= \[Xi] <= Sqrt[1 - 1/M]]
, DirichletCondition[F[u, \[Xi]] == 0.,
umin <= u <= umax && \[Xi] == Sqrt[1 - 1/M]]
};
eqn40 = {(eqn[
u, \[Xi], {\[Mu], \[Eta], S0, \[Tau]c, M,
u0, \[CapitalDelta]u, \[Xi]0, \[CapitalDelta]\[Xi]}] /. {Exp[
x_] -> iExp[x]}) ==
0 + NeumannValue[
0, (\[Xi] == 0 || \[Xi] == 1) && (umin <= u <= umax)] +
NeumannValue[0, (0 <= \[Xi] <= 1) && (u == umin(*||u==umax*))]
, DirichletCondition[F[u, \[Xi]] == 0.,
u == umax && 0 <= \[Xi] <= 1]
};
(* Equation ++Analysis ========================================================*)
Print["getEquations[eqn40]=", getEquations[eqn40]];
(* Reconstruct the parsed differential equation. *)
Print["GetInactivePDE=",
ieqn40LC = GetInactivePDE[getStateObject[eqn40LC]]];
bcData = getBoundaryCondition[eqn40LC];
Print["bcData=", bcData];
(*==============================================================================*)
Print["Solution for Empty Loss Cone"];
fvLC = NDSolveValue[eqn40LC
, F, {u, umin, umax}, {\[Xi], 0, Sqrt[1 - 1/M]}
, "ExtrapolationHandler" -> {Function[0.],
"WarningMessage" -> False}
, Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.002}(*,
"IntegrationOrder"\[Rule]5*)}}
];
Print[fvLC[u, \[Xi]]];
Print[fvLC["ElementMesh"]];
fsLC[u_, \[Xi]_] /; (u < umax && \[Xi] <= Sqrt[1 - 1/M]) =
fvLC[u, \[Xi]];
fsLC[u_, \[Xi]_] /; (u >= umax || \[Xi] > Sqrt[1 - 1/M]) = 0;
Print["Solution for Partially filled Loss Cone"];
fv = NDSolveValue[eqn40
, F, {u, umin, umax}, {\[Xi], 0, 1}
(*,MaxSteps\[Rule]\[Infinity]
,Method\[Rule]{"MethodOfLines",
"SpatialDiscretization"\[Rule]{"TensorProductGrid"}}*)
, "ExtrapolationHandler" -> {Function[0.],
"WarningMessage" -> False}
, Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.002}}}
];
Print[fv[u, \[Xi]]];
Print[fv["ElementMesh"]];
fs[u_, \[Xi]_] /; u < umax = fv[u, \[Xi]];
fs[u_, \[Xi]_] /; u >= umax = 0;
(*===========================================================================*)
(* 3D Plots *)
Print["\n3D plots of Distribution Function"];
Print["\nMagnitudes at {0,\!\(\*SubscriptBox[\(\[Xi]\), \(0\)]\)}=", \
{fvLC[umin, \[Xi]0], fv[umin, \[Xi]0]}];
Print["Maximums=", {(fsLCmax =
NMaximize[{fsLC[u, \[Xi]0], u >= umin && u <= umax},
u]), (fsmax =
NMaximize[{fs[u, \[Xi]0], u >= umin && u <= umax}, u])}];
lbl = StringForm[
"M=``, \!\(\*SubscriptBox[\(v\), \(0\)]\)=``, \
\!\(\*SubscriptBox[\(\[Tau]\), \(c\)]\)=``, \!\(\*SubscriptBox[\(f\), \
\(max\)]\)=``", M, u0, \[Tau]c, fsmax[[1]]];
lblLC =
StringForm[
"M=``, \!\(\*SubscriptBox[\(v\), \(0\)]\)=``, \
\!\(\*SubscriptBox[\(\[Tau]\), \(c\)]\)=``, \!\(\*SubscriptBox[\(f\), \
\(max\)]\)=``", M, u0, 0, fsLCmax[[1]]];
plt3DLC =
Plot3D[Chop@fsLC[Sqrt[vx^2 + vy^2], vx/Sqrt[vx^2 + vy^2]], {vx,
umin, umax}, {vy, umin, umax}
, AxesLabel -> {"\!\(\*SubscriptBox[\(v\), \(z\)]\)",
"\!\(\*SubscriptBox[\(v\), \(\[UpTee]\)]\)", f},
PlotLabel -> lblLC, MaxRecursion -> 4
, MeshFunctions -> {#3 &}, Mesh -> 5
, PlotPoints -> 100
, BaseStyle -> {12}, PlotRange -> pltRange3D, ImageSize -> iSize];
plt3D =
Plot3D[Chop@fs[Sqrt[vx^2 + vy^2], vx/Sqrt[vx^2 + vy^2]], {vx, umin,
umax}, {vy, umin, umax}
, AxesLabel -> {"\!\(\*SubscriptBox[\(v\), \(z\)]\)",
"\!\(\*SubscriptBox[\(v\), \(\[UpTee]\)]\)", f},
PlotLabel -> lbl, MaxRecursion -> 4
, MeshFunctions -> {#3 &}, Mesh -> 5
, PlotPoints -> 100
, BaseStyle -> {12}, PlotRange -> pltRange3D, ImageSize -> iSize];
Print[TableForm[{{plt3DLC, plt3D}},
TableHeadings -> {None, {"Empty LC; reduced diffusion",
"Partially filled LC; reduced diffusion"}}]];
(*===========================================================================*)
(* Density Plots *)
Print["\nDensity Plots of Distribution Function"];
pltDPLC =
Quiet@DensityPlot[
Chop@fsLC[Sqrt[vx^2 + vy^2],
Abs[vx/Sqrt[vx^2 + vy^2]]], {vx, -umax, umax}, {vy, umin, umax}
, FrameLabel -> {"\!\(\*SubscriptBox[\(v\), \(z\)]\)",
"\!\(\*SubscriptBox[\(v\), \(\[UpTee]\)]\)"},
PlotLabel -> lblLC, AspectRatio -> Automatic
, BaseStyle -> {12}, PlotRange -> pltRangeDens,
PlotLegends -> Automatic
, ColorFunction -> "Rainbow"
, MeshFunctions -> {#3 &}, Mesh -> 5
, PlotPoints -> 100
, ImageSize -> iSize];
pltDP =
Quiet@DensityPlot[
Chop@fs[Sqrt[vx^2 + vy^2], Abs@Cos[ArcTan[vx, vy]]], {vx, -umax,
umax}, {vy, umin, umax}
, FrameLabel -> {"\!\(\*SubscriptBox[\(v\), \(z\)]\)",
"\!\(\*SubscriptBox[\(v\), \(\[UpTee]\)]\)"}, PlotLabel -> lbl,
AspectRatio -> Automatic
(*,RegionFunction->Function[{vx,vy,z},Sqrt[vx^2+vy^2]<umax]*)
, BaseStyle -> {12}, PlotRange -> pltRangeDens,
PlotLegends -> Automatic
, ColorFunction -> "Rainbow"
, MeshFunctions -> {#3 &}, Mesh -> 5
, PlotPoints -> 100
, ImageSize -> iSize];
Print[TableForm[{{pltDPLC, pltDP}},
TableHeadings -> {None, {"Empty LC; reduced diffusion",
"Partially filled LC; reduced diffusion"}}]];
(*============================================================================*)
Print["\nDensity Plots of Radial Flux"];
Echo[flux[u, \[Xi], {\[Mu], \[Eta]}], "flux[u,\[Xi]]="];
plt3DfluxLC =
DensityPlot[
flux[Sqrt[vx^2 + vy^2], vx/Sqrt[vx^2 + vy^2], {\[Mu], \[Eta]}][[
1]] /. F -> fsLC, {vx, umin, umax/1}, {vy, umin, umax/1}
, AxesLabel -> {"\!\(\*SubscriptBox[\(v\), \(z\)]\)",
"\!\(\*SubscriptBox[\(v\), \(\[UpTee]\)]\)"}, PlotLabel -> lblLC,
PlotLegends -> Automatic
, ColorFunction -> "Rainbow"
, MeshFunctions -> {#3 &}, Mesh -> 5
, PlotPoints -> 100(*,MaxRecursion->4*)
, BaseStyle -> {12}, ImageSize -> iSize,
PlotRange -> pltRange3DFlux];
plt3Dflux =
DensityPlot[
flux[Sqrt[vx^2 + vy^2], vx/Sqrt[vx^2 + vy^2], {\[Mu], \[Eta]}][[
1]] /. F -> fs, {vx, umin, umax/1}, {vy, umin, umax/1}
, AxesLabel -> {"\!\(\*SubscriptBox[\(v\), \(z\)]\)",
"\!\(\*SubscriptBox[\(v\), \(\[UpTee]\)]\)"}, PlotLabel -> lbl,
PlotLegends -> Automatic
, ColorFunction -> "Rainbow"
, MeshFunctions -> {#3 &}, Mesh -> 5
, PlotPoints -> 100(*,MaxRecursion->4*)
, BaseStyle -> {12}, ImageSize -> iSize,
PlotRange -> pltRange3DFlux];
Print[TableForm[{{plt3DfluxLC, plt3Dflux}}
, TableSpacing -> {3, 3}
, TableHeadings -> {None, {"Empty LC; reduced diffusion",
"Partially filled LC; reduced diffusion"}}]];
(*============================================================================*)
Print["\nDensity Plots of Radial Flux near zero point"];
Echo[flux[u, \[Xi], {\[Mu], \[Eta]}], "flux[u,\[Xi]]="];
plt3DfluxLC =
DensityPlot[
flux[Sqrt[vx^2 + vy^2],
vx/Sqrt[vx^2 + vy^2], {\[Mu], \[Eta]}][[1]] /. F -> fsLC, {vx,
umin, 1}, {vy, umin, 1}
, AxesLabel -> {"\!\(\*SubscriptBox[\(v\), \(z\)]\)",
"\!\(\*SubscriptBox[\(v\), \(\[UpTee]\)]\)"}, PlotLabel -> lblLC,
PlotLegends -> Automatic
, ColorFunction -> "Rainbow"
, MeshFunctions -> {#3 &}, Mesh -> 5
, PlotPoints -> 100(*,MaxRecursion->4*)
, BaseStyle -> {12}, ImageSize -> iSize,
PlotRange -> pltRange3DFlux];
plt3Dflux =
DensityPlot[
flux[Sqrt[vx^2 + vy^2],
vx/Sqrt[vx^2 + vy^2], {\[Mu], \[Eta]}][[1]] /. F -> fs, {vx,
umin, 1}, {vy, umin, 1}
, AxesLabel -> {"\!\(\*SubscriptBox[\(v\), \(z\)]\)",
"\!\(\*SubscriptBox[\(v\), \(\[UpTee]\)]\)"}, PlotLabel -> lbl,
PlotLegends -> Automatic
, ColorFunction -> "Rainbow"
, MeshFunctions -> {#3 &}, Mesh -> 5
, PlotPoints -> 100(*,MaxRecursion->4*)
, BaseStyle -> {12}, ImageSize -> iSize,
PlotRange -> pltRange3DFlux];
Print[TableForm[{{plt3DfluxLC, plt3Dflux}}
, TableSpacing -> {3, 3}
, TableHeadings -> {None, {"Empty LC; reduced diffusion",
"Partially filled LC; reduced diffusion"}}]];
(*============================================================================*)
Print["\nDensity Plots of Flux Norm"];
Echo[Norm@flux[u, \[Xi], {\[Mu], \[Eta]}], "Norm@flux[u,\[Xi]]="];
plt3DfluxNormLC =
DensityPlot[
Norm@flux[Sqrt[vx^2 + vy^2],
vx/Sqrt[vx^2 + vy^2], {\[Mu], \[Eta]}] /. F -> fsLC, {vx, umin,
umax}, {vy, umin, umax}
, AxesLabel -> {"\!\(\*SubscriptBox[\(v\), \(z\)]\)",
"\!\(\*SubscriptBox[\(v\), \(\[UpTee]\)]\)"}, PlotLabel -> lblLC,
PlotLegends -> Automatic
, ColorFunction -> Hue
, MeshFunctions -> {#3 &}, Mesh -> 5
, PlotPoints -> 100(*,MaxRecursion->4*)
, BaseStyle -> {12}, ImageSize -> iSize,
PlotRange -> pltRange3DFlux];
plt3DfluxNorm =
DensityPlot[
Norm@flux[Sqrt[vx^2 + vy^2],
vx/Sqrt[vx^2 + vy^2], {\[Mu], \[Eta]}] /. F -> fs, {vx, umin,
umax}, {vy, umin, umax}
, AxesLabel -> {"\!\(\*SubscriptBox[\(v\), \(z\)]\)",
"\!\(\*SubscriptBox[\(v\), \(\[UpTee]\)]\)"}, PlotLabel -> lbl,
PlotLegends -> Automatic
, ColorFunction -> Hue
, MeshFunctions -> {#3 &}, Mesh -> 5
, PlotPoints -> 100(*,MaxRecursion->4*)
, BaseStyle -> {12}, ImageSize -> iSize,
PlotRange -> pltRange3DFlux];
Print[TableForm[{{plt3DfluxNormLC, plt3DfluxNorm}}
, TableSpacing -> {3, 3}
, TableHeadings -> {None, {"Empty LC; reduced diffusion",
"Partially filled LC; reduced diffusion"}}]];
(*============================================================================*)
Print["\nRadial Flux vs \[Xi] near u=umin"];
us = {1, 2, 4, 10} Max[umin, 0.0001];
plt1DfluxLC =
Plot[(Chop@flux[#, \[Xi], {\[Mu], \[Eta]}][[1]] & /@ us) /.
F -> fsLC // Evaluate, {\[Xi], 0, 1}
, AxesLabel -> {\[Xi], flux}, PlotLabel -> lblLC,
PlotRange -> Automatic
, BaseStyle -> {12}, ImageSize -> iSize, PlotLegends -> us];
plt1Dflux =
Plot[(Chop@flux[#, \[Xi], {\[Mu], \[Eta]}][[1]] & /@ us) /.
F -> fs // Evaluate, {\[Xi], 0, 1}
, AxesLabel -> {\[Xi], flux}, PlotLabel -> lbl,
PlotRange -> Automatic
, BaseStyle -> {12}, ImageSize -> iSize, PlotLegends -> us
];
(*Print[GraphicsRow[{plt3DLC,plt3D}]];*)
Print[TableForm[{{plt1DfluxLC, plt1Dflux}}
, TableSpacing -> {3, 3}
, TableHeadings -> {None, {"Empty LC; reduced diffusion",
"Partially filled LC; reduced diffusion"}}]];
(*====================================================================*)
]
Output
Output by helper functions
Row 1: Plot3D of distribution function
Row 2: DensityPlot of distribution function
Row 3: DensityPlot of radial flux
Row 3-2: DensityPlot of flux norm
Row 4: Plot of radial flux vs. $\xi$ in the vicinity of zero point
FandFFwith exact coefficients as well. Could I give an answer here or it could be better discussed it in a chat? $\endgroup$