dtamayo / reboundx

A library for adding additional forces to the REBOUND N-body integration package
GNU General Public License v3.0
80 stars 60 forks source link

Implementation of Gas Damping Forces #133

Closed PhoebeSandhaus closed 1 month ago

PhoebeSandhaus commented 1 month ago

An implementation of gas damping on particles according to Equation 16 in Dawson et al. 2016. Once you set some depletion factor d which controls the amount of gas present, at each timestep, the code calculates a damping timescale for a particle and applies this timescale to damp both the eccentricity and inclination of the particle.

dtamayo commented 1 month ago

Thanks Phoebe, this looks great! It looks like there are changes to many of the ipynbs...is this because you edited them and added examples to them using your effect, or is that just coming from you running those notebooks locally? If the latter, would you be able to checkout the versions currently on master on this pull request so we don't modify those?

I was also curious whether there's any obvious more specific name we could use for the effect? We now have type I migration, gas dynamical friction, so it seems valuable as we add more to try and be a bit more descriptive in the effect names?

PhoebeSandhaus commented 1 month ago

It was indeed the latter, and I have now checked out the versions from main!

Also, for the name of the function, would "gas_damping_timescale" be better and more descriptive?

dtamayo commented 1 month ago

Hi Phoebe, I'm really sorry for the slow response. This looks great to me--the only remaining question I had was whether you thought that the d parameter should be set on individual particles, or on the effect itself, I.e., if you had two planets, would you want d to be something individual to the two particles, or something that gets set once for the whole disk? Definitely not sure one way or the other--you've thought about this much more than I have! I can merge this right after. Thanks again!

PhoebeSandhaus commented 1 month ago

Hi Dan, no worries! That is a great question!

I've thought about it, and I think allowing the user to hypothetically set the depletion factor for every individual particle would be better so that if someone would maybe want to model gaps/over-densities within a gas disk, they could just give certain particles different values to imitate that effect. But in general, I'd imagine that people would want to just assign the same d parameter to each particle.

Also, I will change my effect name to "gas_damping_timescale" to be more clear about what exactly the function does!

dtamayo commented 1 month ago

Awesome. Also, for the checks to do nothing if e < 1e-7 and i < 1e-7, was that added in response to some numerical or performance issue without it? If not, I would think it might be better to just always do it regardless of the size of the eccentricities or inclinations. I would guess the performance hit is small given that the function needs to be called and the intermediate values calculated, and it would avoid at least one weird discontinuity as parameters change? I'm just adding it to the automated tests and we can merge it

PhoebeSandhaus commented 1 month ago

I can remove those checks since honestly they were mostly to avoid something within my own simulations. And then everything should be good to go!

dtamayo commented 1 month ago

Thanks @PhoebeSandhaus. Did you see my pull request on your fork? I think there's an issue if the user is not using the set of units for the simulation that you're assuming. I made some changes to allow for this, but I also noticed that the resulting damping is significantly slower than you had in your notebook, which I think was due to a units mismatch (units are such a headache!). Could you take a look?

PhoebeSandhaus commented 1 month ago

Okay! I've incorporated your pull request, renamed d_coeff to tau_coeff, and updated the documentation/examples to reflect this change. Let me know if you can think of anything else!

dtamayo commented 1 month ago

Awesome! Were you able to double-check that the smaller damping is correct? (I didn't have a chance to check carefully). I'll add this to the automatic tests in REBOUNDx and once I hear back from you will merge! Thanks again

dtamayo commented 1 month ago

@PhoebeSandhaus I added a test comparing runs with different units, and found a typo in the units I gave above (AU^5/4, not 3/4), so I corrected the notebook and code. I will merge this now to avoid conflicts if we need to make any further changes. Once you confirm if the smaller damping is correct I'll add a new version number and push to PyPI. Thanks!

PhoebeSandhaus commented 4 weeks ago

Hey! I looked at everything and it all seems to be good, except for a couple very tiny things in the Jupyter notebook. Should I make a new pull request or would it just be easier to attach the updated file here and you can manually add it?

dtamayo commented 4 weeks ago

Whatever's easier for you works for me!

dtamayo commented 3 weeks ago

I leave on vacation today, but if you send it to me I'll find a time to add it and update the REBOUNDx version