ascot4fusion / ascot5

ASCOT5 is a high-performance orbit-following code for fusion plasma physics and engineering
https://ascot4fusion.github.io/ascot5/
GNU Lesser General Public License v3.0
29 stars 9 forks source link

Bug in afsi.thermal #109

Closed pbonofig closed 1 week ago

pbonofig commented 4 months ago

I believe there is a small bug in afsi.thermal.

I came across this because I was looking at alphas and neutrons in stellarator geometries where toroidal angle matters.

The default value in afsi.thermal for maxphi = 2pi. However, the phi array is assigned the units of degrees (phi = stuffunyt.deg).

I believe either:

  1. maxphi=360 should be the default or
  2. The phi values are converted from radians to degrees before assigning unyt.deg
pbonofig commented 4 months ago

Actually, I thought this was a units issue, but something else is going on.

When I look at the phi-distribution, it is pretty much at phi=0 and nowhere else, and I am not sure why? I would expect a uniform sampling in phi: Screen Shot 2024-06-03 at 2 54 27 PM Screen Shot 2024-06-03 at 2 54 27 PM

pbonofig commented 4 months ago

Even stranger, when I set minphi=pi/4 or to some other value, the distribution is still only near phi=0.

Also, when I use markergen, to generate markers for input, when nphi=1 for afsi.thermal, then the markers are uniform in phi from [0,360], but the (R,Z)-locations are only from phi=0. When nphi>1, then the phi distribution is like that above (nphi=30).

miekkasarki commented 3 months ago

Could you please clarify if the issue is in the markergen or already at afsi? Neither of those has been properly used for stellarators (at least I haven't) so it's likely that there are bugs related to 3D geometry.

Those distributions are from the afsi output, correct?

pbonofig commented 3 months ago

I think both.

Yes, the above is from AFSI.

When I try maxphi=360, then the plot range goes up to something large like phi~22,000, and there is bit more spread in the above distribution. However, most of the particle source is still near phi=0 like the above.

When I then take this particle source and use it in markergen, I actually get a uniform distribution in phi from [0,360] deg but the (R,Z) locations are almost all sourced at phi=180. I believe it is phi=180 since minphi=0, maxphi=2pi, and nphi=1 sets the edges at [0,360]. It's just a coincidence that my stellarator is two-fold periodic, so the surfaces are identical at phi=0=180. I will also note that this holds true for altered limits.

Adjusting or playing around with minphi, maxphi, and nphi has not produced the expected uniform distributions.

PS - Sorry, I missed the garage today, I forgot it was two back-to-back odd weeks.

miekkasarki commented 3 months ago

I'm able to reproduce this issue with the AFSI tutorial test case. I'll try to fix this within this week.

pbonofig commented 3 months ago

Thanks, please let me know when you have a fix.

I believe the bug similarly applies to afst.beamthermal. I find things are reasonable when one eliminates toroidal dependence (minphi=0, maxphi=2pi, nphi=1) in both the thermal and beam distributions in regards to total source value, but any phi-dependence makes incorrect source profiles.

As I noted in the slack channel, this may propagate to markergen as well, or be an entirely separate bug. Currently, markergen's (R,Z) values do not agree with the phi-angle (at least in 3D geometries).

miekkasarki commented 3 months ago

I think I managed to hunt down the bug. I already pushed it to the develop so could you use that (or cherry-pick commit 27bf73f) to see if the issue is gone? Thanks for spotting this error!

Let me know if there is still a separate issue with the markergen.

pbonofig commented 3 months ago

Thank you! I have just done some quick tests:

  1. My initial tests show much more reasonable toroidal profiles for thermal and beam.thermal afsi sources. I have some lingering questions about the total and thermal sources, though, which I will try to discuss in tomorrow's garage. One such question is: When nphi is increased, so does the total source rate. But what if I want to know the total source rate? How many bins should I choose? I suppose the source rate goes up with nphi because you are simply adding more bins for MC markers?

  2. See the below quick and dirty screenshots. They are of alphas from markergen. The one plot is the corresponding phi per particle id (100k particles). It is clearly uniform in [0,360]. The other plot are the corresponding (R,Z) locations. These only match with the phi=0=180 flux surfaces for my stellarator. The plot range extends to 4.2 in R where I would expect more points at different toroidal angles. So, yes, there also appears to be a bug in markergen regarding the mapping of the (R,Z) locations to the correct toroidal angle.

Screen Shot 2024-06-11 at 12 05 46 PM Screen Shot 2024-06-11 at 12 06 16 PM

miekkasarki commented 3 months ago

Okay, I'll try to have a look at the marker generator as well. The relevant code is rather simple and it's here in case you want to have a look as well.