sirocco-rt / sirocco

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

live or die and extract do not produce the same result #815

Closed kslong closed 1 year ago

kslong commented 3 years ago

Right now, and apparently for some time, the results of live or die and extract differ in both classic and normal modes of operation. The fluxes from extract are systematically lower. The problem appears is readily apparent in agn like models in which continuum emission is dominated by disk emission. The problem does not appear to occur in stellar models, which suggest the issue is with the disk, but the model below has no disk.

A10P0 50_new _new_ld

The parameter file I have been using is below, but it does require James' sed for AGN

new.pf.zip

kslong commented 3 years ago

@jhmatthews Here are some further clues about what is happening:

So it really looks that problem has to do only or at least mainly with photons emerging from the central source. If you think about it, you realize that in the absence of a wind, all of the live or die photons get out, but at angle only half of the the photons do. But if it were solely a geometrical effect then one would expect to see the same thing for a totally low velocity situation as one with high velocities which is not the case. Note that all of the models discussed above were calculated in "classic" mode.

jhmatthews commented 3 years ago

Knox, do we know if this happens with different surface reflection treatments?

kslong commented 3 years ago

No, but the examples I have been running lately, have no wind interactions and so the surface treatment should not matter. We only have one way of treating the angular distribution of emitted photons, which assumes an Eddington approximation.

kslong commented 3 years ago

Following @Higginbottom and @jhmatthews suggestion I compared the created, emitted and 45 deg flux for one of various models I have been looking at xstar, which is a massive BH emitting into a very low mdot wind. For live or die, all of the fluxes agree; for extract, the created and emitted agree, but the extracted flux at 45 degrees is low. So extract is the problem.

jhmatthews commented 3 years ago

I think that's probably right - unless the problem that affects live or die also effects the creation of the Emitted and Created columns (e.g. if it was just to do with packet binning or something).

kslong commented 3 years ago

I have found the error. It does not have anything to do with co-moving frames or special relativity. The problem is that in extract photons were not being reweighted if the origin was PTYPE_AGN . The fixed code is at about line 287 of extract.c and now reads"

For disk and centrol object photons, see Eqn 2.19 Knigge's thesis

 */

  if (itype == PTYPE_STAR || itype == PTYPE_BL || itype == PTYPE_AGN)
  {                             /* It was an unscattered photon from the star */
    stuff_v (pp->x, x);
    renorm (x, 1.);
    zz = fabs (dot (x, xxspec[nspec].lmn));
    pp->w *= zz * (2.0 + 3.0 * zz);
  }

itype==AGN was missing. I suspect these lines of code have not changed for a long time.

Higginbottom commented 3 years ago

Nice one Knox!!

Sent from my iPhone

On 15 Feb 2021, at 22:19, Knox S. Long notifications@github.com wrote:

 I have found the error. It does not have anything to do with co-moving frames or special relativity. The problem is that in extract photons were not being reweighted if the origin was PTYPE_AGN . The fixed code is at about line 287 of extract.c and now reads"

For disk and centrol object photons, see Eqn 2.19 Knigge's thesis

*/

if (itype == PTYPE_STAR || itype == PTYPE_BL || itype == PTYPE_AGN) { / It was an unscattered photon from the star / stuff_v (pp->x, x); renorm (x, 1.); zz = fabs (dot (x, xxspec[nspec].lmn)); pp->w = zz (2.0 + 3.0 * zz); } itype==AGN was missing. I suspect these lines of code have not changed for a long time.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

kslong commented 3 years ago

I have merged the changes regarding the continuum into dev. This does not resolve differences between live or die and extract having to do with anisotropic scattering

kslong commented 3 years ago

Based on a discussion today, @jhmatthews @lazygun37 and I believe the problem with differences between live or die and extract has to do with problems with the interpolation of the velocity gradient tensor and/or dvds max. The current version of Python calculates the both of these values at the corners of cells and interpolates to positions in the cell. The thought is that there is no guaranteed that the values interpolate consistently across a cell. This is causing a problem calculating anisotropic scattering, because we are finding values of dvds that are higher than dvds_max. This in turn leads to incorrect weighting of photons in the extract.

Working with simple models, it is clear that the discrepancy in dvds is substantial sometimes as large as a factor of 2 or more.

As a temporary fix to this problem, and to allow for further testing two changes have been made to Python (currently only the bug815_horns branch). One allows one to avoid calculating the velocity within a cell using the velocity gradient tensor; instead velocities are calculated on the fly. The other allows one to increase dvds_max by a factor that can be set. Both are controlled by parameters in python.h

#define PNORM_FUDGE_FACTOR     5  /*An extra factor used for fudging the velocity factor See #815 */
#define USE_GRADIENTS        TRUE   /*IF true use interpolated velcity gradients to calculate dv_ds */

Setting PNORM_FUDGE_FACTOR to a value other than 1, causes dvds_max to be multiplied by that value in anisowind.c and in extract.c.

Setting USE_GRADIENTS to FALSE means that velocities will be calculated on the fly and the gradient tensor will be avoided.

Both of these choices ameliorate the problem of differences between live or die and extract. Neither of these should be regarded as a final solution, as we need to understand why the interpolation routines are failing as badly as they seem to be.

kslong commented 3 years ago

@jhmatthews With the various changes I have made to Python recently, the changes between live or die and extract have diminished some, but they are still quite apparent in the short AGN model. I had not realized until recently that the biggest problems are at intermediate angles, namely 70 degrees. Furthermore, I had not realized that there were problems both with anisotropic and isotropic scatttering. Here are two figures that show comparisons between live or die and extract, one for anisotropic and one for isotropic scattering (using my latest version of py85d) anisotropic

A70P0 50_test_test_ld

isotropic

A70P0 50_iso_iso_ld

kslong commented 3 years ago

@jhmatthews @lazygun37 The problems that you see in the figures above are very concentrated at an inclination of 70 degrees, which is exactly the inner wind angle inJames' pathological example, which has most of the radiation coming from a small central object of radius about 1e15 with a wind rmax of 1e20.. They go away within at inclination differences of a few degrees. Live or die is averaging over +-2 degrees. This does not necessarily imply that live or die or extract are correct, just that one can understand the differences between the two. The question it begs is whether for this model, do we expect something like this to happen at all, namely to have directions where the the observed flux varied dramatically over a very small angle. (Note that these are all isotropic scattering , classic mode)

Here are a few figures that show how sensitive one is to inclination angle. z.foo.210316.tar.gz

kslong commented 2 years ago

This problem still exists although it may not be as serious as it was. The problem is still concentrated around 70 degrees, as shown in the figure

A70P0 50_disk_ld_disk

It can be reproduced with materials in the attached z.Recreate.220912.tar.gz

kslong commented 2 years ago

I did an experiment in which I extracted only unscattered photons. It's quite clear from this that most of the difference between live or die and extract arises from unscattered photons. Within the noise, the spectra generated with one or more scatters show no real difference. The scattered photons are only 1/10th of the unscattered ones in any event.

Wind.ionization(on.the.spot,ML93,LTE_tr,LTE_te,fixed,matrix_bb,matrix_pow,matrix_est)           matrix_pow
Line_transfer(pure_abs,pure_scat,sing_scat,escape_prob,thermal_trapping,macro_atoms_escape_prob,macro_atoms_thermal_trapping) macro_atoms_escape_prob
Matom_transition_mode(mc_jumps,matrix)             mc_jumps
Atomic_data                                h10_hetop_lohe1_standard80

As shown above this is a macro atom model, but of course we do have photons generated in the wind for the detailed spectra.

(The problem does not change very much is anisotropic scattering is invoked)

kslong commented 2 years ago

The problem is somewhat less apparent but still there in classic mode.

A70P0 50_sc:0_xscat_ld_xscat

(Note that I have generally been assuming that this problem is caused by photons generated in the wind, but of course it could be a pure absorption problem in which for live or die photons it is more likely they will be absorbed along the line of sight.. This is something to check into.)

kslong commented 1 year ago

I have re-investigated where we stand on this issue. The bottom line is that there are small but statistically significant differences (if you run a long enough example) between live or die in both classic and standard mode. It is not obvious where the differences arise. They are not due to the size of the acceptance angle. The larges difference appears at 70 degrees. Simulations done today show a very similar figure to that shown above, except the mysterious spike in the center of the profile is less. This is not a particulary simple way to investigate difference between live or die and extract, and there is no clear path for doing so. The magnitude of the problem is small. So I am closing the issue.