sirocco-rt / sirocco

This is the repository for Sirocco, the radiative transfer code used to winds in AGN and other syatems
GNU General Public License v3.0
27 stars 24 forks source link

Code doesn't deal properly with photons that enter then exit wind in imported models #1025

Closed jhmatthews closed 4 months ago

jhmatthews commented 1 year ago

I'd been seeing strange patterns in photon numbers in an imported JEDSAD model (model files on slack in the mhd-winds channel. The basic issue is that there are big streaks in the photon number: ntot (1) ) I noticed that if you plot the photon positions each call in translate, you get the following plot for a random subset:

traject

This plot shows the locations of the photons just before the call to translate is made. The directions are fine, but some subset of photons passes through the wind correctly (i.e. has lots of points corresponding to passing through each cell) - another subset seems to get to the edge of the wind then bypasses the wind entirely.

I tracked in detail what happens to one photon and this is the result (annotated plot with points from actual code output for photon 30,000:

Screenshot 2023-09-26 at 12 03 27

The basic problem is that, if the photon exits the wind after entering it, the code has no way of allowing it to enter again - because the starting ds it calculates on the call to translate in space after existing the wind is already taking it through the entire wind.

I think this particular problem may only appear in imported models where we don't define wind cones, but I'm not certain. It also begs the question whether.

The fix will clearly have to catch this case and increment ds from some small number to check if you hit the wind, but it's a little confusing to think how best to do this.

jhmatthews commented 1 year ago

One can partly fix this problem by adding a check for when you are between rhomin and rhomax into translate_in_space, and in that case forcing one to find the edge of the wind by gradually incrementing ds rather than using ds_to_wind. This may not the most efficient way to do this, and it also doesn't solve all cases -- here is an updated version of the plot above, where you can see that two of the photons here still show the wrong behaviour:

traject_fix1

jhmatthews commented 1 year ago

The photon below the wind is not interacting with the wind because ds_to_plane returns a negative number and even though the magnitude of this number is less than the distance to the edge of the wind, the code in ds_to_wind ignores the ds if it is less than zero. This is a separate issue which could crop up with other geometries so I'll open a new issue report.

jhmatthews commented 4 months ago

Fixed by #1029