sirocco-rt / sirocco

This is the repository for Sirocco, the radiative transfer code used to model winds in AGN and other systems
https://sirocco-rt.readthedocs.io/en/latest/
GNU General Public License v3.0
30 stars 24 forks source link

Reported radiated luminosity increases linearly with cycles #1005

Closed jhmatthews closed 1 year ago

jhmatthews commented 1 year ago

I couldn't understand why the AGN luminosity was so high in my models - I thought it was an error with upweighting but it's actually just a bug - here's AGN luminosity plotted against cycle as reported in the diag file / output (really just escaping total luminosity from central SED source):

test_lum

I think the error is just that the radiated arrays are not zeroed:

    for (nn = 0; nn < N_ISTAT; nn++)
    {
      z_abs[nn] = 0.0;
      nphot_istat[nn] = 0.0;
    }

    /* loop over the different photon istats to determine where the luminosity went */
    for (nn = 0; nn < NPHOT; nn++)
    {
      z_abs_all += p[nn].w;

      /* we want the istat to be >1 (not P_SCAT or P_INWIND) */
      if (p[nn].istat < N_ISTAT && p[nn].istat > 1)
      {
        z_abs[p[nn].istat] += p[nn].w;
        nphot_istat[p[nn].istat]++;
      }
      if (p[nn].istat == P_ESCAPE)
      {
        radiated[p[nn].origin] += p[nn].w;
      }
      else
        z_else += p[nn].w;
    }

I believe we should be adding a loop over photon orig that zeroes the radiated array before we enter the NPHOT loop.

kslong commented 1 year ago

I really find this bizarre, both because we should have noticed this years ago, and because there is a well defined placed fro rezeroing all the appropriate variables. We need to understand what happened.

jhmatthews commented 1 year ago

I think the reporting was changed slightly - a new radiated array was created around a year ago, so someone just forgot to zero it first, or possibly forgot it isn't re-initialised each cycle because this routine (calculate_ionization) has the loop over cycles inside it.

jhmatthews commented 1 year ago

Should be fixed by 1b972e505f9b85a1a279095