Igor Kotelnikov

I've sent you a message on ResearchGate with my email.

Great! This is what I wanted to ask you about. Keep me informed.

Thid id correct plot.

There is a high probability that one day I will write an article about solving this equation in WM, because the resulting solution refutes some well-established stereotypes in plasma physics. Are you interested in participating in such a publication?

The figure shows the oscillations on a large scale. Their strange shape practically does not depend on the density of the mesh. I have no understanding why they arise.

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.

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.

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.

Q: Is there any other ways to supress the instablity?

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.

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.

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.

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.

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\

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.

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.

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$.

I don't know amplitude of FF at w=0. Above, I took FF[0,[\Xi]]==1. More over I assume that there might be slight variation on $\xi$ even ar w=0.

No, I expect $c e^{-u^2/(2 a^2)}$, that is the soluttion of

Alex, I mean that FF is peaked at $\xi=0$ for $u\lesssim 1$. Here DF approaches isotropic Maxwellian, at least for partially filled Loss Cone.

No, I am insisting, Distribution function is not zero at $\xi=0$, there is not a leak of particles here except for point $u=\xi=0$. Loss cone is located at $\xi>\sqrt{1-1/M}$.

Sorry, Alex. I'm quite sure that the solution is wrong near $v_z=0$, i.e. at $\xi=0$. This is a plane of symmetry; it is entirely located inside the particle confinement region. At a low transverse velocity,$v_{\perp}$ there should be a maximum of the distribution function along the chord $v_{\perp}=const$. Compare my and your DensityPlot. Please, try to restore the Neumann zero condition on the boundaries of the solution area by the angle, at $\xi=0$ and $\xi=1$.

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?

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?

I implemented non-uniform mesh:

@user21 I am trying to invoke ToGradedMesh function introduced in v.13.

@user21. Yes. here are mesh for some of recent runs: NDSolveFEMQuadElement[<3611>] NDSolveFEMQuadElement[<18515>] NDSolveFEMElementMesh[{{4.*10^-6,36.},{0.,1.}},{NDSolveFEMQu‌​adElement[<18515>]}] Perhaps, it is worth to refine mesh locally near u=0? but I don't find how to do that,

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.

@user21: I expect that radial flux near zero velocity, at u -> umin tends to zero. From physics of the problem this would mean that the solution tends to Maxwellian distribution function. I expected that boundary condition NeumannValue[0, (0 <= \[Xi] <= 1) && (u == umin)] provides zeroing of the redial flux. But this is not the case. Best of all it is seen from the right picture in the last row. This picture shows that solution is isotropic in the vicinity of u=umin as expected [we see no dependence on [Xi]]. However it shows thar radial flux does not tend to zero.

Let us continue this discussion in chat.

@AlexTrounev: Done, Note that PlotRange was changed to Automatic.

@AlexTrounev. I removed Chop@ but didn't notice any difference. Then I added "MaxRecursion->4" and again I didn't notice any change. The updated row 3 is now under the old row 3. Please note that radial flux oscillates around 0 near zero of velosity u.