# How do I pose Neumann boundary condition to suppress particles flux into zero point?

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}, $$-\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 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.

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{\frac{2}{\pi }} \sqrt{\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, $$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[iExp][x_] = iExp[x] Derivative[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_] :=
NDSolveProcessEquations[
eq, {u, \[Xi]}, {u, umin, umax}, {\[Xi], 0, Sqrt[1 - 1/M]}
(*,Method\[Rule]{"MethodOfLines",
"SpatialDiscretization"\[Rule]{"FiniteElement"}}*)
][];

getEquations[eq_] := Module[{temp},
temp = getStateObject[eq];
temp = temp["FiniteElementData"];
temp = temp["PDECoefficientData"];
(# -> temp[#]) & /@ {"LoadCoefficients",
"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[]];
lblLC =
StringForm[
"M=, \!$$\*SubscriptBox[\(v$$, $$0$$]\)=, \
\!$$\*SubscriptBox[\(\[Tau]$$, $$c$$]\)=, \!$$\*SubscriptBox[\(f$$, \
$$max$$]\)=", M, u0, 0, fsLCmax[]];

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

(*============================================================================*)
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]}][] /. 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]}][] /. 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]}][] & /@ 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]}][] & /@ 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 4: Plot of radial flux vs. $$\xi$$ in the vicinity of zero point • Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on Mathematica Meta, or in Mathematica Chat. Comments continuing discussion may be removed.
– Kuba
Feb 15 at 10:20
• I have an answer for F and FF with exact coefficients as well. Could I give an answer here or it could be better discussed it in a chat? Feb 15 at 14:31
• @AlexTrounev: bouth ways are fine. Feb 15 at 16:13

I accepted Alexander Trounev's answer and awarded him a bounty in gratitude for his efforts and time spent on solving my question, although it is still far from the final conclusion.

The best result was obtained after replacing the variables $$u\to w=u^2$$. This made it possible to regularize the coefficients in front of all the terms of the equation at zero in energy $$w=0$$, except for the angle diffusion operator $$\partial_\xi\left(\left(1-\xi^2\right)\partial_\xi FF\right)$$. After several attempts to regularize the boundary conditions, as suggested by Alexander in his answer, it was realized that the revealed instability of the solution near zero $$w=0$$ can be suppressed if an additional boundary condition $$\partial_\xi FF=0$$ is imposed at $$w=0$$.

Indeed, the coefficients $$D_{ww}$$, $$A_{w}$$ of the equation at $$FF(w,\xi)$$ and derivatives for $$w$$ are regular at zero and even equal to zero. The coefficients $$D_{\xi\xi}$$ are singular in proportion to $$1/w$$. Apparently, this should mean that $$(1-\xi^2)\partial_{\xi}FF=0$$ must be required for $$w=0$$, otherwise the term that describes the angle scattering has nothing to balance.

Unfortunately, current version of NDSolve' does not allow to set the second boundary condition $$\partial_{\xi}FF[0,\xi]==0$$ for $$w==0$$?

Converting the 2nd-order PDF to a 1st-order PDE system for the vector $$\{FF,Fw=\partial_f, F\xi=\partial_\xi FF\}$$ made it possible to pose all the required BC. However, NDSolve warned that the system is convection-dominated and its solution may be unstable. Indeed, the solution turned out to be too sensitive to the choice of $$wmax$$. Depending on $$wmax$$ it varies from something expected to something wrong.

I still have some ideas to achieve a better solution, but it takes some time to test them. If everything goes as expected, the solution will be posted here.

# Update 2023-03-25

In the end, it turned out that the solution formulated in my question was essentially correct, and the apparent bumpiness in the radial flow plot was the result of incorrect rendering. One way to eliminate this chatter was to add a little bit to the definition of the fs function, which just zero-padded the solution fv of the equation with the NDSolve utility, for example, like this:

fvD1[w_, \[Xi]_] = Derivative[1, 0][fv][w, \[Xi]];
(* The following definitions of fs can be dropped in favour of using \
RegionFunction options in 3D plots *)
fs[w_, \[Xi]_] /; w < wmax = fv[w, \[Xi]];
fs[w_, \[Xi]_] /; w >= wmax = 0;
fsD1[w_, \[Xi]_] /; w < wmax = Derivative[1, 0][fv][w, \[Xi]];
fsD1[w_, \[Xi]_] /; w >= wmax = 0;
Derivative[1, 0][fs][w_, \[Xi]_] = fsD1[w, \[Xi]];

• It is good attempt (+1), but for my opinion it could be better to add some integral norm for F and FF as well, since distribution function should be normalized. In this case we can avoid oscillations on the border. Feb 19 at 3:35
• @AlexTrounev, Could you explain your idea in more details? I managed to set BC FF[0,xi]=FF0=const by removing BC at w=wmax: eqn40w = {(eqnw[ w, [Xi], {[Mu], [Eta], S0, [Tau]c, M, u0, [CapitalDelta]u, [Xi]0, [CapitalDelta][Xi]}]) == 0 + NeumannValue[ 0, ([Xi] == 0 || [Xi] == 1) && (wmin <= w <= wmax)] + NeumannValue[0, (w == wmin)] (,DirichletCondition[FF[w,[Xi]][Equal]0.,w[Equal]wmax]) , DirichletCondition[FF[w, [Xi]] == FF0, w == wmin] }; I expected that this induces some radial ouflux at w=wmax, but it and FF are zero there. Feb 19 at 6:44

As it well known the final form of numerical solution for this problem is M.F=B, where matrix M should be inverted. In a case of homogeneous Neuman boundary conditions only we have singular matrix and regularization problem discussed here. To avoid singularity we can update NeumannValue[0,...] as in a case of the Tikhonov regularization. This method works good for equation with the variable w=u^2, while for the equation with u we have some improvement as well but not final solution. The code for F is given by

Needs["NDSolveFEM"]

ClearAll[iExp]
iExp[x_] /; (x >= -676.) = Exp[x]
iExp[x_] /; (x < -676.) = 0
Derivative[iExp][x_] = iExp[x]
Derivative[iExp][x_] = iExp[x]

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

inactiveFei[
u_, \[Xi]_, {\[Mu]_, \[Eta]_}] = -Dei .
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

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[Pi/4]
, \[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_] := NDSolveProcessEquations[ eq, {u, \[Xi]}, {u, umin, umax}, {\[Xi], 0, Sqrt[1 - 1/M]} (*,Method\[Rule]{"MethodOfLines", "SpatialDiscretization"\[Rule]{"FiniteElement"}}*) ][]; 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]}) == NeumannValue[( F[u, \[Xi]])/umin^2, u == umin && (0 <= \[Xi] <= Sqrt[1 - 1/M])] + NeumannValue[ F[u, \[Xi]]/umin^2, (\[Xi] == 0) && (umin <= u <= umax)] , DirichletCondition[F[u, \[Xi]] == 0., u == umax && 0 <= \[Xi] <= Sqrt[1 - 1/M]] , DirichletCondition[ F[u, \[Xi]] == 0., \[Xi] == Sqrt[1 - 1/M] && umin <= u <= umax] }; eqn40 = {(eqn[ u, \[Xi], {\[Mu], \[Eta], S0, \[Tau]c, M, u0, \[CapitalDelta]u, \[Xi]0, \[CapitalDelta]\[Xi]}] /. {Exp[ x_] -> iExp[x]}) == NeumannValue[( F[u, \[Xi]])/umin^2, u == umin && (0 <= \[Xi] <= 1)] + NeumannValue[ F[u, \[Xi]]/ umin^2, (\[Xi] == 0 || \[Xi] == 1) && (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.001}(*, "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.001}}} ]; 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[]]; lblLC = StringForm[ "M=, \!$$\*SubscriptBox[\(v$$, $$0$$]\)=, \ \!$$\*SubscriptBox[\(\[Tau]$$, $$c$$]\)=, \!$$\*SubscriptBox[\(f$$, \ $$max$$]\)=", M, u0, 0, fsLCmax[]]; 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->{"Subscript[v, z]","Subscript[v, \ \[UpTee]]"},PlotLabel->lblLC,PlotLegends->Automatic ,ColorFunction\[Rule]Hue ,MeshFunctions->{#3&},Mesh\[Rule]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->{"Subscript[v, z]","Subscript[v, \ \[UpTee]]"},PlotLabel->lbl,PlotLegends->Automatic ,ColorFunction\[Rule]Hue ,MeshFunctions->{#3&},Mesh\[Rule]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]}][] & /@ 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]}][] & /@ 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"}}]]; Print["\nRadial Flux vs u at \[Xi]=\[Xi]0"]; plt1DfluxLC = Plot[(Chop@flux[u, \[Xi]0, {\[Mu], \[Eta]}][]) /. F -> fsLC // Evaluate, {u, umin, umax} , AxesLabel -> {u, flux}, PlotLabel -> lblLC, PlotRange -> Automatic , BaseStyle -> {12}, ImageSize -> iSize]; plt1Dflux = Plot[(Chop@flux[u, \[Xi]0, {\[Mu], \[Eta]}][]) /. F -> fs // Evaluate, {u, umin, umax} , AxesLabel -> {u, flux}, PlotLabel -> lbl, PlotRange -> Automatic , BaseStyle -> {12}, ImageSize -> iSize ]; Print[TableForm[{{plt1DfluxLC, plt1Dflux}} , TableSpacing -> {3, 3} , TableHeadings -> {None, {"Empty LC; reduced diffusion", "Partially filled LC; reduced diffusion"}}]]; (*====================================================================*) ] In a case of FF form with w=u^2 we update boundary conditions as follows eqn40wLC = {(eqnw[ w, \[Xi], {\[Mu], \[Eta], S0, \[Infinity], M, u0, \[CapitalDelta]u, \[Xi]0, \[CapitalDelta]\[Xi]}] /. {Exp[ x_] -> iExp[x]}) == NeumannValue[ FF[w, \[Xi]]/ umin, (\[Xi] == 0 || \[Xi] == 1) && (wmin <= w <= wmax)] + NeumannValue[ FF[w, \[Xi]]/ umin, (0 <= \[Xi] <= Sqrt[1 - 1/M]) && (w == wmin)] , DirichletCondition[FF[w, \[Xi]] == 0., w == wmax && 0 <= \[Xi] <= Sqrt[1 - 1/M]] , DirichletCondition[FF[w, \[Xi]] == 0., wmin <= w <= wmax && \[Xi] == Sqrt[1 - 1/M]] }; eqn40w = {(eqnw[ w, \[Xi], {\[Mu], \[Eta], S0, \[Tau]c, M, u0, \[CapitalDelta]u, \[Xi]0, \[CapitalDelta]\[Xi]}] /. {Exp[ x_] -> iExp[x]}) == NeumannValue[ FF[w, \[Xi]]/ umin, (\[Xi] == 0 || \[Xi] == 1) && (wmin <= w <= wmax)] + NeumannValue[ FF[w, \[Xi]]/umin, (0 <= \[Xi] <= 1) && (w == wmin)] , DirichletCondition[FF[w, \[Xi]] == 0., w == wmax && 0 <= \[Xi] <= 1] };  • Thanks, Alex. I have a few days to review your solution before the bounty expires. In the past, I had to deal with Tikhonov regularization and SVD. Since then, I have been left with the feeling that regularization is more of an art than a science. Could you explain the choice of regularizing functions? Have you tried other options? Why, for example, did you modify the BC and not the diffusion coefficients? Do you think your method will work when passing to a non-stationary problem? Why did you regularize the BC for chi=0 and chi=1? After all, this didn't seem to be a problem, did it? Feb 16 at 8:15 • First of all, this regularization works well for FF with exact coefficient only. For the model with F it is a good attempt to the final solution. Second, this is really art to find out the right function for regularization Mathematica FEM matrix. But nevertheless it works. Concerning nonstationary case we need to test my approach. Other options are solutions with using LDG method and wavelets method. Do you like these methods or FEM is ultimate? Feb 16 at 8:58 • What confuses me most about your solution is that the regularization led to a decrease in the maximum of the distribution function by hundreds of times compared to the solution without regularization. Similarly, I was previously confused by that the addition of the second derivative with respect to v caused this maximum to be increased hundreds of times. This is critically important for the stability criterion with respect to DСLС, one of the two main kinetic instabilities in an anisotropic plasma. Maybe we will continue chatting, otherwise I'm afraid that we will be sent there forcibly? Feb 16 at 9:21 • But in your model there is no any flux on the border$\xi =0$and only very sharp function Source$ with maximum at u0 = 5, \[Xi]0 = Cos[45. Degree]`. How you suppose to get something on the border $\xi =0$? Feb 16 at 12:29
• Do you mean Maxwell distribution $c u^2 e^{-u^2/(2 a^2)}$? Actually it comes to zero at $u\rightarrow 0$. Feb 16 at 14:21