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

Very optically thick CV model failing (due to population inversions in first ionization cycle) #1019

Closed kslong closed 1 year ago

kslong commented 1 year ago

Yasuke has a model that is failing at the end of the first ionization cycle. The model has a very large amount of scattering.

No. of photons which have scattered n times.     The max number of scatters seen was 1999
4.33e+06 1.96e+06 9.96e+05 6.23e+05 4.35e+05 3.23e+05 2.54e+05 2.07e+05 1.74e+05 1.5e+05

Most of the scatters arise apparently from electron scattering, as


Number of photons resonantly scattering n times.  The max number of scatters seen was 23
9.51e+06 1.85e+06 6.37e+05 2.61e+05 1.19e+05 5.78e+04 2.78e+04 1.46e+04 8.08e+03 4.85e+03
3.04e+03 2.02e+03 1.45e+03 971      667      457      334      226      149      125
91       57       34       33

Huge numbers of photons are hiitting the disk

No of photons and their fates
!!PhotFate: 0        0        8.03e+06 3.23e+05 4.2e+04  0        1.47e+05 3.95e+06 0        0

Photons contributing to the various spectra
                      Inwind    Scat      Esc       Star      >nscat    err       Absorb    Disk     Sec        Adiab(matom)
             Created  0         0         8.03e+06  3.23e+05  4.2e+04   0         1.47e+05  3.95e+06  0         0
            WCreated  0         0         0         0         0         0         0         0         0         0
             Emitted  0         0         8.03e+06  0         0         0         0         0         0         0
              CenSrc  0         0         7.62e+03  0         0         0         0         0         0         0
                Disk  0         0         4.32e+06  0         0         0         0         0         0         0
                Wind  0         0         3.7e+06   0         0         0         0         0         0         0
             HitSurf  0         0         0         3.23e+05  0         0         0         3.95e+06  0         0
           Scattered  0         0         3.99e+06  0         0         0         0         0         0         0

All this leads to a population inversions, and the program exiting

!!python: Number of ionizing photons 1.71379e+64 lum of ionizing photons 4.22039e+53
Error: normalise_simple_estimators: Huge w 3.39e+19 in cell 30 trad   2.75e+04 j 3.51e+32
init_freebound: Creating recombination coefficients
init_freebound: Creating recombination emissivities between 0.000000e+00 and 1.000000e+50 in structure element 0
Error: sobolev: VERY BAD population inversion in cell 30: d1 1.19556e-05 d2 3.58669e-05 g1 1 g2  3 freq 2.51286e+15 f 0.0476 frac_upper 0.78881
Error: sobolev: VERY BAD population inversion in cell 30: d1 1.19556e-05 d2 3.58669e-05 g1 1 g2  3 freq 2.58816e+15 f 0.00557 frac_upper 0.78881
Error: sobolev: VERY BAD population inversion in cell 30: d1 1.19556e-05 d2 3.58669e-05 g1 1 g2  3 freq 2.67325e+15 f 0.00135 frac_upper 0.78881
Error: sobolev: VERY BAD population inversion in cell 30: d1 1.19556e-05 d2 3.58669e-05 g1 1 g2  3 freq 2.6817e+15 f 0.000987 frac_upper 0.78881
Error: sobolev: VERY BAD population inversion in cell 30: d1 1.19556e-05 d2 3.58669e-05 g1 1 g2  3 freq 2.68961e+15 f 0.003 frac_upper 0.78881
Error: sobolev: VERY BAD population inversion in cell 30: d1 4.11128e-05 d2 6.16692e-05 g1 4 g2  6 freq 2.64139e+15 f 0.0116 frac_upper 0.198842
Error: sobolev: VERY BAD population inversion in cell 30: d1 4.11128e-05 d2 4.11128e-05 g1 4 g2  4 freq 3.10786e+15 f 0.0105 frac_upper 0.132561
Error: sobolev: VERY BAD population inversion in cell 30: d1 4.11128e-05 d2 4.11128e-05 g1 4 g2  4 freq 3.14362e+15 f 0.0231 frac_upper 0.132561
Error: sobolev: VERY BAD population inversion in cell 30: d1 377.35 d2 1132.05 g1 1 g2  3 freq 2.76563e+15 f 0.109 frac_upper 1.22315
Error: sobolev: tau is >1e-3 and nu < gu/gl * nl. Exiting.
Para: --------------------------------------------------------------------------
Aborting rank 0001: exiting all processes with error 1

The .sig file indicates that it is possible this happening in the wind_update stage

Fri Sep 15 10:43:25 2023   4141.1 NOK                  Photon transport completed
Fri Sep 15 10:43:35 2023   4150.5 NOK                  Start wind update0

The pf file is below:

05030.pf.txt

It is not immediately obvious what is leading to the situtaion, allthough there is quite a bit of heading of the disk, which has been set to absorbing.

kslong commented 1 year ago

This problem seems to arise because at the end of the first cycle w (the ratio of the observed intensity to that expected for a BB is so large that we have too much stimulated emission.)

In the example, above the problem goes away if the mass loss rate is slightly lower, or if the starting T is 1e5, instead of 4e4.

I believe this to be a transient effect due to a "poor" starting model. I have added comments at the point where the program exits for the issue on how to get around the problem. The comments suggest raising the temperature of the input model, when the problem arises, and if not to reopen this issue

I am closing this issue.