« first day    last day (20 days later) » 

Igor Kotelnikov
Igor Kotelnikov
02:11
Let's summarize the (interim) results.

1. The best results of solving the equation are obtained after the transition from the spherical coordinate system {u,\theta,\alpha} to the Cartesian {w=u^2, \xi=\cos\theta} (the solution does not depend on \alpha).

2. In these variables, `NDSolve` gives the correct solution for the empty loss cone model. At the same time, my initial boundary conditions were used: FF=0 at w=wmax, FF=0 at \xi=sqrt{1-1/M} (at the boundary of the loss cone) and the null Neumann condition Flux_w=0 at w=-wmin, and you can even take wmin=0. The solution found in this way h…
3. In the model of a partially filled cone of losses in such Cartesian coordinates, we see that the radial flux of particles Flux_w = D_{ww} \partial_w FF +A_w FF to zero in energy w monotonically decreases from the maximum at the ejection energy w=wmax to almost zero at w=wmin. This means that the boundary conditions are set correctly, namely: FF=0 at w=wmax, and zero Newman at w=wmin, $\xi=0$, $\xi=1$.
4. According to the type of solution graphs (that is, according to the FF distribution function), the solution found meets all expectations, but suspicion was aroused by the fact that it has a giant peak near zero in energy. No known publication has reported such a peak, however, there are almost no publications with a model of a partially filled cone of losses.
5. However, on the graph of the radial component of the Flux_{ww} flow, an oscillation was observed in the area of approximately $u\lesssim 0.5$. It indicates the instability of the solution near zero energy.
6. The use of a smaller grid only slightly reduced the width of the unstable zone and also slightly increased the wavelength of the oscillations. At the maximum, I used 800,000 nodes of a rectangular uneven grid.
getMesh[{wMin_, wMax_, wCount_,
wMinimal_}, {\[Xi]Min_, \[Xi]Max_, \[Xi]Count_}] :=
Module[{mesh, mesh, mesh},
mesh =
ToGradedMesh[
Line[{{wMin}, {wMax}}], <|"Alignment" -> "Left",
"ElementCount" -> count, "MinimalDistance" -> minimal(*,
"GradingRatio"\[Rule]1.1*)|>];
meshY =
ToGradedMesh[
Line[{{\[Xi]Min}, {\[Xi]Max}}], <|"Alignment" -> "Uniform",
"ElementCount" -> \[Xi]Count(*,"GradingRatio"\[Rule]1.5*)|>];

mesh = ElementMeshRegionProduct[meshX, meshY]
]
8. To generate an uneven grid, the getMesh function was used, which is defined above. The following command
mesh = getMesh[{wmin, wmax, 1000, 1./100}, {0, 1., 50}];
generated a grid with 50,000 nodes.
Igor Kotelnikov
Igor Kotelnikov
03:04
Computation and drawing of all plots took less than 5 minutes. Next 2 pictures show radial flux profile along the angle of particles injection. We see no oscillations of the plot for the model of empty loss-cone and see such oscillations of the same plot for the model of partially filled LC\
user image
user image
Igor Kotelnikov
Igor Kotelnikov
03:19
9. Reducing the lifetime of the particle $\tau_c$ in the loss cone by 10 times, from $\tau_c=0.01$ to $\tau_c=0.001$ practically did not change the graph of the radial projection of the flow, but reduced the peak amplitude on the distribution function $f_{max}$ from 9 to 3.4, which is already enough close to the value of 0.65 in the empty loss cone model.
user image
Igor Kotelnikov
Igor Kotelnikov
03:46
10. Sorry, I just have noted that in fact the above pictures were drawn for the case where zero NeumannValue BC at w=wmin was subsitituted with the BC proposed by Alex Trounev, i.e.
eqn40w = {(eqnw[
w, \[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) && (wmin <= w <= wmax)] +
NeumannValue[FF[w, \[Xi]]/umin, (0 <= \[Xi] <= 1) && (w == wmin)]
, DirichletCondition[FF[w, \[Xi]] == 0.,
w == wmax && 0 <= \[Xi] <= 1]
}
11. Initially, Alexander Trutnev proposed to regularize the solution by replacing all zero Neumann boundary conditions with non-zero ones, for example, NeumannValue[FF[w, [Xi]]/umin, (0 <= [Xi] <= 1) && (w==wmin)]. In this way, it was possible to destroy the oscillator on the Flux_{ww} chart. However, the distribution function has acquired incorrect features. In particular, in the limit w->0, it did not approach isotropic.
Igor Kotelnikov
Igor Kotelnikov
04:08
12. Replacing the Neumann boundary condition at all boundaries where this condition was used made it possible to suppress the oscillation. However, replacing the boundary condition at $\xi=0$ led to an incorrect physical result, as mentioned above. On the other hand, replacing the boundary condition only at w=wmin did not solve the problem of suppressing instability.
Q: Is there any other ways to supress the instablity?
Alex Trounev
Alex Trounev
04:44
Igor, can you formulated your problem correctly as possible? You discussed solution properties like isotropic distribution around zero while problem looks as incorrect. It could be better to put isotropic distribution as boundary condition and add integral equation to fix norm of solution (number of particles is fixed).
 
3 hours later…
Igor Kotelnikov
Igor Kotelnikov
07:22
Alex, Alex, let's think together. The coefficients $D_{ww}$, $A_{w}$ of the equation at FF 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$ for $w=0$, otherwise the term that describes the angle scattering has nothing to balance.
Is it possible to set the second boundary condition $\partial_{\xi}FF[min,\xi]==0$ for $w==wmin$? Maybe it will help somehow? I'll think about the integral equation.
 
10 hours later…
Igor Kotelnikov
Igor Kotelnikov
17:26
Converting the 2nd-order PDE to a 1st-order PDE system has yielded encouraging results. But it's too late. I'm going to bed.

« first day    last day (20 days later) »