Runaway temperature with matrix_est and macro atom line transfer

[Higginbottom]

In certain, possibly extreme, situations, the electron temperature can be driven to probably unphysically high temperatures. It is not clear exactly why this happens but it requires matrix_est as the ionization scheme, and macro-atom line transfer.

This is a thin shell model that reproduces the 'problem'

input.pf.txt

It produces a thin shell, with density of n_H=1e11cm^-3, and an ionization parameter of U=1.1e8 xi=1.6e10.

With various different parameter, we ge the following behaviour of temperature with cycles...

the 1e5 etc represents the number of photons..

matrix_pow gives a reasonable and very stable temperature

matrix_est is OK for standard line transfer, but a little noisy. Interestingly, increasing the number of photons does not change the noise...

matrix_est with macro atom line transfer shows a runaway, increasing the number of photons makes no odds. The runaway is caused by large PI heating, which does not appear to be reflected in a suitably increased ionization state, so more heating can take place until Compton cooling eventually comes to the rescue...

[kslong]

211107 - ksl - Verified that this problem still exists and is straightforward to reproduce.

[kslong]

@Higginbottom @jhmatthews The reason this error is occurring is that the in macro atom mode, the photoionization rates are all 0. This is because the values of xplasma->ioniz[mm] are all 0, where the photonization rates are calculated in matrix_ion.c, near line 140. The reason xplasma->ioniz[mm] are all 0 is that these are incremented in the subroutine radiation, which is called from translate in wind, near line 500, but if you look at the code there, you will see that radiation is not called in macro_atom mode.

   if (geo.rt_mode == RT_MODE_MACRO)
   {
     /* In the macro-method, b-f and other continuum processes do not reduce the photon
        weight, but are treated as as scattering processes.  Therefore most of what was in
        subroutine radiation for the simple atom case can be avoided.  
      */
     if (geo.ioniz_or_extract == CYCLE_IONIZ)
     {
       /* Provide inputs to bf_estimators in the local frame  */
       stuff_phot (p, &phot_mid);
       move_phot (&phot_mid, 0.5 * ds_current);
       observer_to_local_frame (&phot_mid, &phot_mid_cmf);
       ds_cmf = observer_to_local_frame_ds (&phot_mid, ds_current);
       bf_estimators_increment (&w[p->grid], &phot_mid_cmf, ds_cmf);
     }
   }
   else
   {
     radiation (p, ds_current);
   }

What do you all think we should do about this.

[kslong]

This problem still exists.

[jhmatthews]

Fixing this properly is a little fiddly and would require adding the radiation code to bf_estimators_increment. And not a priority given how rarely it is used - plus, as we move towards more macro-atoms this mode will become even less important within this line transfer mode. However, the code should produce an error and quit if you try to run with matrix_est and macro-atoms, so I should put this error catch in and then close this issue.