Photons are hitting the disk and the star and being lost when they should be reflected

[kslong]

When using the reflection option, photons should not be lost when they hit the star or disk, but this is happening in TDE and presumably other models.

The issue is related to #584, and at least in the case of photons hitting the star to a recent change made in walls (to deal with photons which travel long distances) and were not being caught as hitting the star.

Part of the problem is that walls is being used for a variety of things. It presumes one wants to find the boundary of something (the disk, or star mainly) given two photons, one starting with where you are and one a proposed new location. If you hit something along the way, it tells you and reduces the distance the old photon has travelled. But this only works if the photon has not been redirected in the mean time.

[kslong]

I made a change, that modifies how walls checks for a star. If the photon direction has changed between old and new, I assume we only care if the new position is inside or outside the star. I regard this as a kluge at present, but it does reduce but not totally eliminate phtons that are recorded as hitting the disk.

I have begun to track down what happens in the remaining cases of this. It looks like the problem is a repositioning one. similar to the one we had for the disk previously. A photon that was in_wind, is being repositioned near line 591 of trans_phot. And it is being repositioned into the star. It looks like there is no attempt in this case to extract ourselves from this problem (as there is for a disk). My guess at this point is that repositioning for the disk is also not completely effective, but I have not verified this.

[kslong]

As noted above, as the code exists in dev, the check in walls, for whether a photon has intersected the wind looks like this:

  s = ds_to_sphere (geo.rstar, pold);

  if ((r = dot (p->x, p->x)) < geo.rstar_sq)
  {
    /* Then the photon is inside the star */
    hit_star = TRUE;
  }

  else if (s < VERY_BIG && ds_to_sphere (geo.rstar, p) == VERY_BIG)
  {
    /* then we hit the star somewhere in between */
    hit_star = TRUE;
  }

  if (hit_star == TRUE)
  {
    stuff_phot (pold, p);
    move_phot (p, s);
    stuff_v (p->x, normal);
    return (p->istat = P_HIT_STAR);
  }

The second if statement was added to eliminate a problem #635 having to do with getting different fluxes at different inclination angles in a spherical wind. It was intended to deal with the case where a photon has gotten from an old position on one side of the star to a new position on the other side of the star. But this assumes that the old photon and new photon are on the same path. But that is not the case, when a resonant scatter occurs as the new photon has a new direction.
To correct this, I modified the code above the change the if statement to read

   if ((r = dot (p->x, p->x)) < geo.rstar_sq)
  {
    /* Then the photon is inside the star */
    hit_star = TRUE;
  }

  else if (s < VERY_BIG && ds_to_sphere (geo.rstar, p) == VERY_BIG && dot (p->lmn, pold->lmn) > 0.99)
  {
    /* then we hit the star somewhere in between */
    hit_star = TRUE;
  }

This greatly reduced the energy lost to hitting the star. Following on this, I made changes to assure that when walls was called from trans_phot_single, that the pold was updated before calling walls. This completely eliminated the problem of photons hitting the star and the disk, but ...

[kslong]

I then checked to see what the spectra looked like in the Ed's small TDE model. This is the spectrum for a restarted version of tde_cv_slightly_smaller.

And this is what my updated version produces

[kslong]

I have compared the old and new versions of the code. The changes I made to try to fix the problem affect the SNe model in the regression tests, but not in a very significant way the other examples that are there. This may be what one would suspect because these models have far more limited backscatters onto the star or disk than does the tde model, which stresses everything that happens close to the disk and star.

[saultyevil]

So unfortunately it doesn't look this issue is fully fixed just yet. I'm still seeing problems with photons hitting the disc and photons being lost to DFUGDE.

!!python: Total photon luminosity before transphot 7.326901838412e+43

 Ion. Cycle 17/20 of tde_cv : Photon          0 of     402728 or    0.0 per cent 
 Ion. Cycle 17/20 of tde_cv : Photon      40272 of     402728 or   10.0 per cent 
 Ion. Cycle 17/20 of tde_cv : Photon      80544 of     402728 or   20.0 per cent 
 Ion. Cycle 17/20 of tde_cv : Photon     120816 of     402728 or   30.0 per cent 
 Ion. Cycle 17/20 of tde_cv : Photon     161088 of     402728 or   40.0 per cent 
 Ion. Cycle 17/20 of tde_cv : Photon     201360 of     402728 or   50.0 per cent 
 Ion. Cycle 17/20 of tde_cv : Photon     241632 of     402728 or   60.0 per cent 
 Ion. Cycle 17/20 of tde_cv : Photon     281904 of     402728 or   70.0 per cent 
 Ion. Cycle 17/20 of tde_cv : Photon     322176 of     402728 or   80.0 per cent 
Status of    348968 changed from 0 ti 7 on reposition
 Ion. Cycle 17/20 of tde_cv : Photon     362448 of     402728 or   90.0 per cent 
 Ion. Cycle 17/20 of tde_cv : Photon     402720 of     402728 or  100.0 per cent 

!!python: photon transport completed in: 7465.015581 seconds
Error: trans_phot: 3 photons were lost due to DFUDGE (=5.0000e+07) pushing them outside of the wind after scatter
!!python: Total photon luminosity after transphot  7.326901838412e+43 (absorbed/lost  0.000000000000e+00). Radiated luminosity 7.303054501654e+43
!!python: luminosity lost by adiabatic kpkt destruction 0.000000000000e+00 number of packets 0
!!python: luminosity lost to low-frequency free-free    0.000000000000e+00 number of packets 0
!!python: luminosity lost by being completely absorbed  0.000000000000e+00 
!!python: luminosity lost by too many scatters          2.384733673836e+41 
!!python: luminosity lost by hitting the star           0.000000000000e+00 
!!python: luminosity lost by hitting the disk           2.699144213100e+32 
!!python: luminosity lost by errors                     0.000000000000e+00 
!!python: luminosity lost by the unknown                0.000000000000e+00 
spectrum_create: Fraction of photons lost:  0.00 wi/ freq. low, 0.00 w/freq hi

This was created for cycle 17/20, for the following parameter file tde_cv.pf.txt with filling factor f = 0.01. This also happened for the same model, but with a filling factor f = 0.1.

Obviously this has made a huge difference though. I don't see any photons hitting the star, and we aren't losing ~10% of luminosity to the star anymore. It looks like it's just the occasional photon hitting the disc.

[saultyevil]

Also, it does look like this change affects the results. These following spectra are for filling factors f = 0.1 with solar abundances,

AGN (i = 75, 85 degrees)

CV (i = 60, 85 degrees)

Looks like we might need to tweak our best bet models a little bit and double check the code is doing what we want now!

(sorry for the poor quality of the plots, I'll post higher resolution ones onto Evernote)

[lazygun37]

Do we understand why this is yet? What's the optical depth in these models along various sightlines?

Is the convergence much better now? And is luminosity conserved?

C

On 25/11/2019 13:16, Edward Parkinson wrote:

Also, it does look like this change affects the results. These following spectra are for filling factors f = 0.1 with solar abundances,

AGN (i = 75, 85 degrees)

[Screenshot from 2019-11-25 13-13-49]https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fuser-images.githubusercontent.com%2F19627279%2F69543453-8bdffa00-0f85-11ea-91f9-cddf52f742c3.png&data=01%7C01%7CC.Knigge%40soton.ac.uk%7C1a4b407049744ab3cbba08d771a9a9f3%7C4a5378f929f44d3ebe89669d03ada9d8%7C0&sdata=DdjTjLrlcbNgtD36fcBPmU52dEeq2XtqYvQ3oA3v%2FMw%3D&reserved=0

CV (i = 60, 85 degrees)

[Screenshot from 2019-11-25 13-12-56]https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fuser-images.githubusercontent.com%2F19627279%2F69543441-8682af80-0f85-11ea-8dbc-d6ea57dbce89.png&data=01%7C01%7CC.Knigge%40soton.ac.uk%7C1a4b407049744ab3cbba08d771a9a9f3%7C4a5378f929f44d3ebe89669d03ada9d8%7C0&sdata=O2SHH7Ml4Q8QQYIdUSdxnJdP8gCzjPfkbiHRXswc63U%3D&reserved=0

Looks like we might need to tweak our best bet models a little bit and double check the code is doing what we want now!

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fagnwinds%2Fpython%2Fissues%2F645%3Femail_source%3Dnotifications%26email_token%3DABEXZAFYMAAOACD4JF7LIKTQVPF2JA5CNFSM4JPZXPN2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEFCLGTI%23issuecomment-558150477&data=01%7C01%7CC.Knigge%40soton.ac.uk%7C1a4b407049744ab3cbba08d771a9a9f3%7C4a5378f929f44d3ebe89669d03ada9d8%7C0&sdata=wbqd5Rr76VrdMzCTI8QuHXeBHX7Iy45AoSMd9LuwrTI%3D&reserved=0, or unsubscribehttps://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FABEXZAAACIMUP36D72JYRYTQVPF2JANCNFSM4JPZXPNQ&data=01%7C01%7CC.Knigge%40soton.ac.uk%7C1a4b407049744ab3cbba08d771a9a9f3%7C4a5378f929f44d3ebe89669d03ada9d8%7C0&sdata=31UcwQOrLN%2Fq62xrO0YuXb%2F99mDyQh9To6x1S1W8Bzc%3D&reserved=0.

--

Professor Christian Knigge Physics & Astronomy University of Southampton Southampton SO17 1BJ

[saultyevil]

So, for the CV model at least...

Convergence is better with the proposed fix. Sight lines through the wind do not seem to be as optically thick. But the optical depth along the disc plane has increased.

The luminosity is conserved now it seems. What goes in prior to trans_phot comes out after trans_phot apparently (see above Python output in my previous comments). Though I'm slightly confused by the fact that we are losing 1e41 ergs to too many scatters, but Python is reporting absorbed/lost = 0.000000000000e+00. Is this something we expect for macro atoms? I should also note that it does look like the number of electron scatters has increased. Previously, on average per MPI processes ~204 photons for 833300 photons in total would be killed due to too many scatters. But now ~83582 photons for 833300 total are killed. The number of resonant scatters has seemed to decrease as well.

I still need to check things out for the AGN and spherical models. But straight away I can at least say that convergence is better.

[saultyevil]

Okay - I think I may have an idea why the results are so different...

In this commit, the call to the scattering function was commented out and in subsequent commits it was removed from the code. So my guess is that photons are never actually having their directions changed.

To test if this is happening, I ran a quick reduced version of the equatorial TDE model and dumped out the each photon to file, and plotted a few of the photon trajectories. They looked something like this,

Basically, photons never changed direction and either reached MAXSCAT or escaped from the system with their emission direction. So, next steps are to put the scattering function back in and see if we get out lines back :-).

[kslong]

Please be more specific about where in the code this happens.

Sent from my iPad

On Nov 26, 2019, at 5:46 AM, Edward Parkinson notifications@github.com wrote:



Okay - I think I may have an idea why the results are so different...

In this commithttps://github.com/agnwinds/python/commit/6ce2d974cb4ec60dd045700d1913e7c603781fe2#diff-86b7ee4c3c3441647eb814a4f99df768, the call to the scattering function was commented out and in subsequent commits it was removed from the code. So my guess is that photons are never actually having their directions changed.

To test if this is happening, I ran a quick reduced version of the equatorial TDE model and dumped out the each photon to file, and plotted a few of the photon trajectories. They looked something like this,

[photon_2]https://user-images.githubusercontent.com/19627279/69627211-230a8780-1042-11ea-995e-29eae5de1a35.png

Basically, photons never changed direction and either reached MAXSCAT or escaped from the system with their emission direction. So, next steps are to put the scattering function back in and see if we get out lines back :-).

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/agnwinds/python/issues/645?email_source=notifications&email_token=AATJ4VPFOCMKMQJ5YFY2DEDQVUEBTA5CNFSM4JPZXPN2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEFFXJBI#issuecomment-558593157, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AATJ4VIKF5CVGOOWQF54VY3QVUEBTANCNFSM4JPZXPNQ.

[saultyevil]

Sorry - maybe the next links should make it clearer.

You should be able to see in this commit that scatter was commented out, and then in this commit was actually deleted from the source code.

[saultyevil]

This function call was around line 467 in trans_phot.c in the trans_phot_single function.

[saultyevil]

The spectra from the log_spec_tot file after the first ionisation cycle looks more promising, but I can't be sure just yet. This is for the AGN style TDE model by the way, so won't be great to compare with the spectra Knox previously posted.

Note that some photons are still hitting the surface. It looks like the error condition for issue #584 is being triggered (I need to double check this is working as intended now) and curiously some photons are still defiantly hitting the star. What's confusing is for an example photon I pulled out which hits the star, nres = -1. So it is not being pushed into the star because of reposition. It's not clear to me what is missing now...

[kslong]

That is definitely a mistake on my part, and should be restored. I thought I was just eliminating errors that never came up. But I will be surprised if that solves the problem as this was supposed just to be cleanup.

[saultyevil]

Okay so, adding the scatter function looks like it bought back the absorption lines. Phew! Below is a synthetic spectra for an unconverged AGN style model - I only ran with a small number of photons for a few cycles.

Where previously a similar model looked like this,

But this means that the original problem has unfortunately not been fixed. There are now a lot of errors reporting,

Error: walls: the direction p and pold differ

Which I can understand. This pretty much happens because a photon can scatter (here) and have its direction changed. Then in the subsequent call to walls () after reposition () (here), p and pp are naturally going to have different directions which I believe is triggering the error condition,

if (dot (p->lmn, pold->lmn) < 0.99)
{
  Log ("Error: walls: the direction p and pold differ\n");
}

I think this is related to the original problem, specifically this line,

else if (s < VERY_BIG && ds_to_sphere (geo.rstar, p) == VERY_BIG)

What I think is happening is that pold was pointing towards the star, hence s = ds_to_sphere (geo.rstar, pold) returns a value s < VERY_BIG. But, p may have a new direction which is not pointing towards the star anymore, so ds_to_sphere (geo.rstar, p) == VERY_BIG. I think I've understood that right? If the photon doesn't intersect the sphere on it's current path, then VERY_BIG is returned?

Knox said in an earlier post that this was added to fix issue #635 for when a photon passes straight through a star (i.e. because of translate_in_wind I guess?) when it shouldn't. However, it seems to me that with this logic, it's possible for a photon which was pointing towards the star, but is then scattered and now points away from the star to be falsely identified as hitting the star. If my thought process is a bit of a mess and unclear, I have no problem making a crude diagram to better show what I think is happening.

CK and I thought a little harder about this and we think we can simplify this a bit. Since we know where the photon was before trans_phot with pold, we can calculate the distance from pold to hitting the star using ds_to_sphere (). We also previously calculated in calculate_ds () the distance that photon is going to move, ds. So, we thought it would be possible to implement a check like this instead,

if (dot (p->x, p->x) < geo.rstar_sq || p->ds > s)

That way this first checks to see if the photon is in the star and then checks to see if the photon moved further than s, which is the distance from the old position to the surface of the star. We think this should work and catch cases where a photon passes straight through the star - but of course would need testing with Mandy's parameter files from issue #635. This means we don't need to worry about the directions of the old and new photon. What do you think @kslong? Can you think of any reason why this wouldn't work?

[kslong]

I am glad that this restored the lines. I was too quick to try to get rid of lines that I thought had not effect.

I may not have done this correctly, but what was intended, was that by the time you got to walls, one had the old and new photon going in the same direction. The reason that photons have a different direction before they get to walls, unless you do something is that they are redirected by the scatter that occurs earlier.

I am not sure your change to the condition for hitting the sphere will work. What the problem was that on extraction, you had photons that originated on one side of the star say from (-1,0,-1) . In the direction of (1,0,1) and these were getting through. I would like you to verify that you are getting the same result with the old and new version of the condition.

Knox

[saultyevil]

I ran the extract parameter files from #635 using the version of the code in my fork, and I believe there is no excessive flux problem. I've re-created Mandy's plots and it seems reminiscent of the final plots she was happy with,

I still need to run the regression tests to see how this effects our standard models, as well as test the clumpy TDE models again - and a general documentation effort as well.

[kslong]

I still would like to see a straight up comparison between the two approaches

[saultyevil]

Okay, so I've done the same tests for a couple of commits now and overplotted my approach and the various previous approaches. There are lots of plots below!

I tested my approach against the original approach using commit 10da17e1aaf11dfc8381104552ed30fd4d14e1fa. This should be the latest commit with scatter added back in and with the original approach.

Here are the same parameter files, but now using a commit before all of the changes related to this issue, 85cfefc0dc4e4feaf665540d4bfd75d169878724.

Finally, here are the same parameter files but with a comparison before the fix. It's not super clear, but it's instructive enough to show that the issue is there. I think this shows that the problem can be handled by the new approach, especially if you look carefully at the plot with the smaller wavelength range.

[kslong]

OK, if you are convinced they are all the same, and that Mandy’s problem remains solved, I am happy, and we should assume your modification has worked. If we are all happy then we should merge the changes into dev.

K

[saultyevil]

So it looks like photons hitting the disc comes back to the reposition problem of Issue #584. It looks like there are multiple odd scenarios where a photon can be pushed underneath the disc and trying to catch all of these different cases was difficult.

I think there are a few solutions,

• A complicated if statement and possibly an extra flag for when a photon has been repositioned.
• Change to what we pass to walls after reposition. For example, we would probably need a photon pp_before_reposition to boundary check, but then this still requires additional code to deal with bad repositions.
• Catch the bad smax < 0 case, but allow P_HIT_DISK, P_HIT_STAR and put the reflection code into a function which can be called after reposition if a photon hit the disc or star.
• Catch the bad cases and decide not to reposition these photons which cause trouble.

Option 2 does still require a solution on how to deal with the bad photons. For me, I think options 3 and 4 are both reasonable. I have actually implemented option 4 because this was the easiest. Nothing hits the disc or star now and regression tests seem ok - the changes are small, apart from the 1d_sn model which never seems to have small changes now. Though, option 4 would require a decent chunk of testing to see if we are ever in the regime where the Sobolev interaction region that we would accidentally interact with the same line again. I suspect not for most models? But I guess this might eventually cause problems where there are lots of resonance interactions by the disc plane, for example if try an even denser TDE model.

[kslong]

Please place a .pf file here that shows the problem.

[saultyevil]

The original .pf I used to test this problem was using CNO abundances - which I haven't been able to push to the atomic data branch yet. However, it appears to also happen with solar abundances. I have confirmed that we see this problem on the latest version of dev. You will note that the error is something like,

Status of 30015 changed from 0 to 7 after reposition

The number of these errors is small as most errors are from from 0 to 12, which is the istat for reposition error. Either way, here is the pf.

tde_cv_disk.pf.txt

[kslong]

OK, I confirm this is still happening, but as Ed says, rarely.

[kslong]

@saultyevil In the context of trying to deal with vertically extended disks #700, I may have found the cause of this problem. walls assumes that the two photons it are given, the old position and the new position are moving in the same direction. However on resonant scattering walls was being called when the old position did not have the same direction as the new photon (scattered) photon. I don't believe this affects a flat disk, but it certainly affected hits to the central object.

I will add that the logic of trans_phot on a scatter is not very transparent, presumably as result of various attempts to fix problems. We should try to figure what is really happening and make a cleaner fix if possible. Or if we have a clean fix, we should delete the tests.

[kslong]

Even after the changes that have been made to vertically extended disks, this error has not disappeared. It looks like one can get the error both to 3 (HIT_STAR) and 7 (HIT_DISK)

[kslong]

I actually did not understand this. tde_cv_disk.pf has a flat disk, and so the problem is unlikely to be related to a vertically extended disk. This is a separate problem from #700

[kslong]

This should be closed with recent updates to reposition and in the way walls works.