Add photoevaporation modules#4
Conversation
… old timescale method which will work for high constant mass loss rates and may be useful for testing. Also has the first attempt at weighting for the optically thin regime. Jupiter Mass added to constants for convenience. Photoevaporative timescale added to driver to avoid situations where more mass than exists can be lost in the timestep. Extra features on disc designed to account for half empty cells and ambient UV field.
… not yet removing entrained material
…0. Implementation of this, along with using integrated mass. Tracking additional disc properties
…entrainment working. Work to tackle problems of over-rapid radial drift and diverging dust mass fractions. Currently setting the minimum size to zero again, which avoids the drift problem, but we still add a check that it is drifting slower than the ffree-fall velocity or terminal velocity. Attempt to model pressure gradient commented out. To solve the diverging densities, we reset the dust mass fraction to 1 should we exceed it, though there is a commented out attempt to diagnose the pproblem and do better.
…d ffmpeg into the python script. Updated run_many to iterate over several, use nice and run in the background, waiting for batches to finish before moving onto the next.
…ct of this on the drift timescale and nan errors. Changed file saving to logarithmic spacing in line with the plotting. Considered whether to entrain dust when truncating - commented out for now
… scripts and save more regularly. Estimating depletion times and fate fractions.
…so that other things can call it more easily.
…photoevaporating area is a hole.
… done within, where appropriate by calling functions in other modules to return a value.
|
I think I've made the structural changes you suggested. |
|
Thanks for this, it looks great. I need to make some changes to get the FRIED grid to work with my install scripts (which have changed in the last couple of years). Once I've done that, I'll let you know |
|
Sounds great – glad to hear it’s all coming together!
|
|
Do you have time to chat today about the viscous velocity implementation? |
|
It's ok, I think its fine |
|
Ok, I am free any time between now and 5 if anything else does come up or you want to check the implementation anyway.
Andrew
|
|
Can you try installing my latest version and checking that your internal / external PE results look sensible? |
|
The only siginificant changes I've made are:
|
|
Yes I can try this now
|
|
On the second of these – the prescription of Dipierro+2018 used 1/rho_g to define their v_visc, so that it corresponds to the velocity of the gas in the limit that there is no coupling - this is why I passed the gas density. But I see you’ve multiplied by 1+eps in the dust module to account for this which is probably a simpler implementation – I wonder if this should be 1/1-eps though?
|
|
I think multiplication is correct. The torques (forces) are the same but the total mass moved is more, so I get a lower velocity. The 1+eps should correct that |
|
I agree that the lower velocity results - my issue is that their v_visc is independent of the normalisation of rho_g so at fixed total rho, as rho_d -> rho_tot, v_visc shouldn’t change, but the actual radial velocity of the gas (their eq 20) does decrease.
From that equation taking St=0 and vp=0 (so everything moves at one speed determined by viscosity as per your argument) and accounting for their definition of epsilon being different: 1/(1+rho_d/rho_g) = rho_g/rho_tot = (1-eps) -> ur = (1-eps) * v_visc. Hence dividing by (1-eps) (which for smaller eps ~ * (1+eps) ) should also make the velocity lower?
I’ll have to run and check I get sensible results with changing this in the cavity/gap when eps->1.
|
|
I think all my tests look sensible.
With regards to the question of the 1+eps factor, I ran the same X-ray disc with my original way of calculating the viscous term, and the change you made. Mostly I don’t think it really gets to the point where eps~1 and these differences would matter except in the gap where there is very little dust anyway and no gas to provide drag/couple to. As the final clearing of the outer disc goes on, eps does start getting larger. Since the dust is trapped anyway it seems to have little effect, but the feedback on the gas is visible – with your implementation the growth of the hole radius is lightly slower which keeps the mass loss rates a bit higher.
These differences might be more important in the case of the debris discs Cathie’s student is going to be working on, so I think it is important that we agree which is right…
Just to let you know I identified a couple of bugs - one which was sometimes stopping the code correctly identifying a hole and another where I hadn’t changed when I reset the hole radii if one that later turns out to not be the correct hole disappears (it sometimes gets confused by the outer disc depending on the interplay between the Sigma and Sigmadot profiles).
I’ll update these on my version and put in a pull request to fix them in your master.
|
|
I see my mistake now, and agree with you that 1+eps should be replaced with 1 / (1-eps). Do you want to make that change. Thanks for fixing bugs, too. |
|
Ok I'll make that change along with the others. Do you suggest I make the denominator 1+1e-300-eps for stability?
|
|
It might help, was the small number (1e-300) needed in your original implementation? |
|
You could just use eps_g as already defined in the code instead of 1-eps.sum() |
|
I didn’t need it explicitly because I was passing the gas density which presumably is similarly pinned to a small number.
But yes, this is the most sensible idea I think. I’ll make the change now and test it.
|
|
I’m having some issues with this at the moment – when the gas density is sufficiently low in the gap (say eps_g ~ 10^-17) then when the code calculates eps_av, it loses this precision that eps is not quite 1, and return eps_av=1 -> average eps_g=1e-300. This means that the correction factor is a couple of hundred orders of magnitude too high, which leads to insanely high velocities, grinding everything to a halt.
I think if I calculate SigG_av and eps_g explicitly rather than from the dust it ought to eliminate this problem, so I’m going to test that out now.
|
|
That should work since its essentially what you were doing before. |
|
Ok that worked – I’ve put in a new pull request with the edits. I think it’s good we got this right but actually the differences I highlighted before seemed to be arising in the hole identification because it was using the Picogna prescription to clear up anything left in the gap which had an outer edge so lead to a non-contiguous photoevaporating regions – I’m not sure why the bug only showed up when we changed the implementation of the gas velocities but it should all be fixed now.
From: Richard Booth <notifications@github.com>
Sent: 10 November 2020 18:38
To: rbooth200/DiscEvolution <DiscEvolution@noreply.github.com>
Cc: Andrew Sellek <ads79@cam.ac.uk>; Author <author@noreply.github.com>
Subject: Re: [rbooth200/DiscEvolution] Add photoevaporation modules (#4)
That should work since its essentially what you were doing before.
—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub<#4 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AJXA3ZYMKHXIJ4THWFU44MLSPGCAZANCNFSM4OPZYQ6Q>.
|
2 participants
Additional external and internal photoevaporation modules and edits to existing modules to work with these. Upgrade of the dust evolution to use the prescription from Dipierro et al. 18 for the dust-gas relative velocity.