Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Photons absorbed in macro atom model of a dense CV wind #1024

Open
kslong opened this issue Sep 24, 2023 · 13 comments
Open

Photons absorbed in macro atom model of a dense CV wind #1024

kslong opened this issue Sep 24, 2023 · 13 comments

Comments

@kslong
Copy link
Collaborator

kslong commented Sep 24, 2023

Photons are not normally supposed to be absorbed in macro atom models, and yet in a dense CV wind that is exactly what is happening. These photons are being reported as absorbed and in 87d they appear carry a small amount of energy, but in fact the value that is being reported is not the original energy but the the final energy carried by the photons that is being reported.

foo.pf.txt

@kslong
Copy link
Collaborator Author

kslong commented Sep 24, 2023

This is part of a print out of some variables atter an ionization cycle.

The lines labelled stat refer to status of photons after the transport, e.g escaped 2, absorbed = 6. What is reported is the number of photons, the final weight of those photons, and the initial weight of the photons.

The lines labelled rad, refer to emitted photons. but here the numbers are associated with the origin (and whether they are macro atom photons are not.

The main point of this is to show that the absorbed photons actually represent a large fraction of the total luminosity of the system.

!!python: photon transport completed in: 98.738525 seconds
xxx Sun Sep 24 14:59:25 2023    106.6 NOK                  Photon transport completed
XXX stat        0            0         0.000e+00       0.000e+00
XXX stat        1            0         0.000e+00       0.000e+00
XXX stat        2       282086         2.247e+35       2.753e+35
XXX stat        3            0         0.000e+00       0.000e+00
XXX stat        4            0         0.000e+00       0.000e+00
XXX stat        5            0         0.000e+00       0.000e+00
XXX stat        6       134580         6.060e+28       1.637e+35
XXX stat        7            0         0.000e+00       0.000e+00
XXX stat        8            0         0.000e+00       0.000e+00
XXX stat        9            0         0.000e+00       0.000e+00
XXX stat       10            0         0.000e+00       0.000e+00
XXX stat       11            0         0.000e+00       0.000e+00
XXX stat       12            0         0.000e+00       0.000e+00
XXX rad        0        4.278e+32       4.284e+32
XXX rad        1        0.000e+00       0.000e+00
XXX rad        2        1.781e+35       1.785e+35
XXX rad        3        0.000e+00       0.000e+00
XXX rad        4        0.000e+00       0.000e+00
XXX rad        5        0.000e+00       0.000e+00
XXX rad        6        0.000e+00       0.000e+00
XXX rad        7        0.000e+00       0.000e+00
XXX rad        8        0.000e+00       0.000e+00
XXX rad        9        0.000e+00       0.000e+00
XXX rad       10        1.142e+32       2.199e+32
XXX rad       11        0.000e+00       0.000e+00
XXX rad       12        4.602e+34       9.611e+34
XXX rad       13        0.000e+00       0.000e+00
XXX rad       14        0.000e+00       0.000e+00
XXX rad       15        0.000e+00       0.000e+00
XXX rad       16        0.000e+00       0.000e+00
XXX rad       17        0.000e+00       0.000e+00
XXX rad       18        0.000e+00       0.000e+00
XXX rad       19        0.000e+00       0.000e+00

@kslong
Copy link
Collaborator Author

kslong commented Sep 25, 2023

@jhmatthews @lazygun37

The problem is occurring here

Hardware watchpoint 9: pp.w

Old value = 3.7634096718185621e+26
New value = 7.4169601326339181e+25
scatter (p=0x7fffffffd750, nres=0x7fffffffd6f4, nnscat=0x7fffffffd6f0) at resonate.c:1102
1102	          macro_gov (p, nres, 2, &which_out);   //routine to deal with kpkt
2: prob_kpkt = 0.19708085963040745
3: p->w_orig/p->w = 166.94810355415294
4: *p = {x = {-2962297871.4855661, 1874136098.8262272, 647767978.13935328}, lmn = {0.39856551291945708, -0.64088882677881231, -0.65605414686695429}, freq = 16545003481571256, freq_orig = 968201694504242.75, w = 7.4169601326339181e+25, w_orig = 1.2382474282799913e+28, tau = 0, istat = P_SCAT, frame = F_LOCAL, nscat = 6, nrscat = 0, nmacro = 5, nres = 200132, line_res = -3, nnscat = 1, grid = 60, origin = PTYPE_STAR_MATOM, origin_orig = PTYPE_STAR, np = 1, path = 2940407395.672895, ds = 57694.811732947303}
(gdb) where
#0  scatter (p=0x7fffffffd750, nres=0x7fffffffd6f4, nnscat=0x7fffffffd6f0) at resonate.c:1102
#1  0x000055555557b100 in trans_phot_single (w=0x55559f896630, p=0x7fff872e70a8, iextract=0) at trans_phot.c:430
#2  0x000055555557a56c in trans_phot (w=0x55559f896630, p=0x7fff872e7010, iextract=0) at trans_phot.c:131
#3  0x00005555555eb4a6 in calculate_ionization (restart_stat=0) at run.c:225
#4  0x0000555555559b70 in main (argc=2, argv=0x7fffffffe368) at python.c:653

The problem can be avoided by turning off "upweighting", which here is downweighting photons regularly.

@kslong
Copy link
Collaborator Author

kslong commented Sep 25, 2023

The "upweigthing" for bf is carried out a few lines earlier in the code.

else if (phot_top[*nres - NLINES - 1].macro_info == FALSE || geo.macro_simple == TRUE)
      {
        /* Simple ion case.  It's a bf interaction in a calculation involving macro-atoms, 
           but this photoionization x-section is is associated with one of the levels
           a simple atom, not one that is part of a full blown macro atom.

           (Alternatively we are treating all atoms in a simplfied mode)

           Need to make decision about making a k-packet. Get the fraction of the energy
           that goes into the electron rather than being stored as ionisation energy: this
           fraction gives the selection probability for creating a k packet. It's given by 
           (1 - ionization edge frequecy / photon frequency)

           Note:  In some cases, one can obtain a negative prob_kpkt below.  This happens because
           we use the  mean frequency along the path to calculate continuum opacities
           in calculate_ds, whereas here we have the exact frequency at a particular point.
           This leads to a situation where the co-moving frequency can be less than
           the edge frequency.  See issues #436 and #867.

           If the error is large, it suggests that the length of that a photon is allowed
           to travel in a single step is too large.

         */

        prob_kpkt = 1. - (phot_top[*nres - NLINES - 1].freq[0] / freq_comoving);

        if (prob_kpkt < 0)
        {
          /* only report an error for a negative prob_kpkt if it's large-ish in magnitude. see #436 discussion */
          if (prob_kpkt < -1e-2)
          {
            Error ("scatter: kpkt probability (%8.4e) < 0 for phot_top %d, zeroing\n", prob_kpkt, *nres - NLINES - 1);
            Log ("scatter: photon edge frequency: %8.4e, comoving (observer) frequency %8.4e %8.4e\n", phot_top[*nres - NLINES - 1].freq[0],
                 freq_comoving, p_orig.freq);
          }
          prob_kpkt = 0.0;
        }

As written it always reduces the photon weight. How is the energy supposed to come out @jhmatthews

@jhmatthews
Copy link
Collaborator

See the docs here: https://agnwinds.readthedocs.io/en/dev/physics/macroatoms.html?highlight=upweight#bound-free-continua-of-simple-atoms - there's another place where we upweight.

Note also that we record and report the amount of energy being extracted/deposited in the bound-free simple ion pool. I think the two bits of code are in kpkt() and scatter() -- although one of the comments is a bit confusing there and seems to have got mixed up.

@kslong
Copy link
Collaborator Author

kslong commented Sep 26, 2023

@jhmatthews The situation described in Readthedocs is not what is happening. In this case, the problem is that the photons are being downweighted so much from multiple scatterings that they eventually count as absorbed. Note that I have tried Christian's suggestion for running a model to convergence (75%) without upweighting, and then running another cycle with it on. While we don't loose 50% of the flux as we do when you run with upweighting on, we do lose 25%, which I think should be concerning (even if we make not upweighting the default.)

Is there anywhere that we have a full description of how this is all supposed to work; the readthedocs page really does not explain this anywhere. If this does not exist, it would be extremely useful if one were created.

@kslong kslong changed the title Photons absorbed in macro atom model of a dense CV sind Photons absorbed in macro atom model of a dense CV wind Sep 26, 2023
@jhmatthews
Copy link
Collaborator

I posted the Readthedocs page because you asked above how the energy is supposed to come out, not because it explains this issue which has only just come up in this particular form. The cautionary note about numerical problems is indeed a different issue (although related).

In terms of how all this is supposed to work, the description on RtD, whilst doubtless flawed, is my best attempt to explain that. Perhaps there are some things missing but I can't think what exactly.

@kslong
Copy link
Collaborator Author

kslong commented Sep 29, 2023

There was a suggestion made by @ssim that this problem might be due to the fact that we do not use the Milne relation to calculate recombination rates might result in the problem we are seing, because the recombination rates might be higher than predicted. In that case, one might eliminate the problm if we used the Milne relation instead. @jhmatthews suggested we could test this by removing the separate recombination rates from the program. So I modified the data file we have been using h20_hetop_standard80.dat, removing these lines
''''
data/atomic/recomb.dat
data/atomic/chianti_dr.dat
data/atomic/chianti_rr.dat
data/atomic/Badnell_GS_RR.dat
data/atomic/di_dere.dat
''''
but this had no effect. We are still losing a lot of energy into the ion pool

The logs say that the Milne relation is being used, so presumably this is a good test of this hypothesis.

@kslong
Copy link
Collaborator Author

kslong commented Sep 30, 2023

We document the flow of energy in the log files for the entire wind. For the original model and the one which is supposed to follow the Milne relation this is what one gets:

!! report_bf_simple_ionpool: Total flow into: 1.0706e+35 and out of: 3.0572e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.0153e+35 and out of: 9.7539e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9328e+35 and out of: 5.5755e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9107e+35 and out of: 3.7056e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9225e+35 and out of: 3.0034e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9375e+35 and out of: 2.5415e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9348e+35 and out of: 2.3325e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9404e+35 and out of: 2.1314e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9451e+35 and out of: 2.0195e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9428e+35 and out of: 1.9931e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9571e+35 and out of: 1.8773e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9604e+35 and out of: 1.9484e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9579e+35 and out of: 2.0849e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9628e+35 and out of: 2.0314e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9689e+35 and out of: 2.1613e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9715e+35 and out of: 2.0468e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9681e+35 and out of: 2.1715e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9849e+35 and out of: 2.1706e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9851e+35 and out of: 2.1656e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9806e+35 and out of: 2.5205e+34 bf_simple ion pool

In this case the energy in is typically 10 x the energy out. Note there is one line for each of the 20 cycles

BUT if one runs a macro atom model, without any macro atoms you get something quite different

!! report_bf_simple_ionpool: Total flow into: 1.7142e+36 and out of: 1.6677e+36 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 9.2296e+35 and out of: 7.4214e+35 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.2051e+36 and out of: 1.0200e+36 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 4.9550e+38 and out of: 4.9610e+38 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.8208e+40 and out of: 3.1459e+40 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.0955e+41 and out of: 5.4736e+40 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.9381e+39 and out of: 3.1784e+39 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 6.0569e+37 and out of: 6.5355e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 5.5339e+36 and out of: 6.0763e+36 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 3.6363e+38 and out of: 3.6920e+38 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.7500e+37 and out of: 3.2789e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.0501e+37 and out of: 1.1638e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 6.9131e+37 and out of: 6.7843e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.1363e+37 and out of: 3.3468e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.5875e+36 and out of: 2.9579e+36 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.3165e+36 and out of: 2.2220e+36 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 7.4902e+36 and out of: 7.3923e+36 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.2931e+37 and out of: 2.2221e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.3686e+37 and out of: 2.3999e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 3.1149e+37 and out of: 3.1298e+37 bf_simple ion pool

The totals are higher now, and there is variability. But in this case the energy in is equal to the energy out.

@jhmatthews @ssim Is this expected? Is it possible that something is wrong when a simple atom is photoionized but a macro atom is selected for recombination.

@kslong
Copy link
Collaborator Author

kslong commented Oct 3, 2023

@jhmatthews On branch bug1024_lost, I added more tracking of what bf process are doing. The file shows (I think) what fb process have absorption and emission. The first column is the fb no, the second is the elment involved, the third is the ion state and the 4th and 5 columns are the number of times one absorbed by that fb process and the number of decays by that process.

Some of it makes sense, we see lots of re-combinations to H and He. What does not make sense is the huge number of interactions with NV and OVI, since both are almost absent in the wind.

foo.txt

@jhmatthews
Copy link
Collaborator

I agree that's confusing. I guess it's possible that there are some sizable minority of cells where there is a decent amount of NV and OVI, since I think this file is a cumulative one across the whole wind? It's not clear to me exactly what we mean here by "number of times one is absorbed by fb process" i.e. whether it means r-packet activations of a simple-atom via that bound-free process? If so, the chance of it happening is set just by the bound-free opacity, I think.

The post before about macro-atoms is less surprising to me, because removing the macro-atom data is changing how the H and He continua are treated, so I might well expect the problem to increase (the downweighting can also take place in the H and He continua if they are treated as simple ions, rather than being done "self consistently" within the macro-atom).

@kslong
Copy link
Collaborator Author

kslong commented Oct 3, 2023

The numbers here are, as James guessed the total for the wind, but as far as I can tell the OVI and NV concentrations are very very small. If one does a diff between the bug1024_lost branch, it is easy to see where I did the counting.

I guess my question about the difference between results with and without macro atoms is for the H and He macro atom case, do these create a situation where necessarily we are creating an imbalance. Photons that undergo bound free absorption necessarily lose the ionization energy, but when they come back through a macro atom recombination, they may not get it back. I know the upweighting which occurs early is expected to compensate for that, but it does not seem to.

@kslong
Copy link
Collaborator Author

kslong commented Oct 3, 2023

Well, my previous comment was incorrect. Indeed OVI and NV are the dominant versions of these species, as the disk in this model is very hot 10**5 K. So it rather looks like the following is happening. BF transitions are dominated by ionization of ions like NV and OVI, but they are recombining mainly into H and He. Photons are losing ionization energy when they are absorbed, but the boosting that occurs is insufficient to make up for the energy that they lose.

The question that needs answering it seems to me now is why the upweighting does not, on average at least, balance this all out.

@ssim
Copy link
Collaborator

ssim commented Oct 4, 2023 via email

kslong added a commit that referenced this issue Oct 16, 2023
* Provides additional diagnostics of energy flow into and out of the ion pool (see #1024)
* Version changed to 87f due to the fact that there are additional variables in the plasma structure
* There should be no actual changes to the way photons are processed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

3 participants