Ode

[tigerchenlu98]

The up-to-date code for self-consistent spin, tidal & dynamical EoMs

[dtamayo]

I've actually been working on putting in these transformation functions in too, though I think it might make sense to put them in REBOUND. I think everything is working, but I am finishing the documentation and an example. I'll put in a REBOUND pull request tomorrow that we can all look at and decide what the best way to organize it might be!

[tigerchenlu98]

Great to hear that the transformation functions are working! Went ahead and removed my attempt to do it - I agree REBOUND is a more natural place for it!

[tigerchenlu98]

@dtamayo @hannorein - all examples (python & C) have been updated with the new syntax!

I'd like to fully run through some of the longer C examples to double check they still give the expected results - I can't quite get to that tonight, but I don't expect anything to have gone wrong.

One suggestion - I do think it would be useful to directly have helper functions that:

• align the simulation to the invariant plane
• rotate the spin vector into the planet frame since right now both of those are 3-liners in the C code, if I'm not mistaken. I personally ended up needing both of those specific transformations quite a bit, so I think others would get good use out of them too. Happy to work on these functions if you think they'd be helpful!

[dtamayo]

Hi Tiger, this looks great! I really like the Earth-Moon example you put together.

@hannorein and I talked about that quite a bit, and I think he convinced me. It's better to have a general API that can be used for anything you could want to do, rather than a really large code base with many different functions that the user has to go through the documentation to understand and get the syntax for. I think it's also expected for things to be a bit more verbose in C. Certainly everything can still be an intuitive one-liner in Python.

I think we're really close to merging!

• This is actually the perfect effect to profile REBOUNDx with, because there are so many force evaluations and getting and setting of parameters. I followed your analysis, and reproduced your result that most of the time is being spent in rebx_get_type, because it needs to do a strcmp with every param in the list of registered parameters.

Thankfully, that only over needs to happen when we create a new parameter, which is only done initially and isn't a big performance hit. Unfortunately, it was needlessly being done every timestep. I think I never noticed before because there weren't that many get/set params, and our list of registered params started out much smaller--so this would have got worse over time as more params get added. Thanks so much for catching that! After the refactor, the code is 2x faster, the same as with your set_spin_param (so I've now removed that).

It's an interesting question whether we should go beyond that and get rid of the linked lists. I think just over half the time on these kinds of runs is spent in set and get param, so it seems right around the uncomfortable boundary where we might want to do it! Hanno expressed some interest in possibly looking into it. I can't focus on it short-term, so let's not worry about it for this release.

• I'm not sure rebx_spin_potential is complete--shouldn't there be a contribution to the potential from the spins? I think this is likely a carryover oversight from the previous CTL implementation, but I also modified it to only require the source to have a non-zero mass (OK if it doesn't have sigma, r or k2), and the target to have non-zero mass, radius and k2. Those are the only things we use in calculate_spin_potential.

• I wanted to make sure I'm not getting a convention wrong, but the ipynbs calculated

sx = s sin(obl) sin(phi) sy = s sin(obl) cos(phi)

rebound.spherical_to_xyz is doing the opposite:

sx = s sin(obl) cos(phi) sy = s sin(obl) sin(phi)

I thought phi would always be measured from the x axis. Was that a typo or did I get that wrong?

[hannorein]

Oh wow, good catch with the 2x speedup. And yes, I'm happy to look into speeding up these linked list operations. I think there are a couple of easy things one could try. Basically, just caching various pointers. Given how rarely the cache will get invalidated (when parameters/particles get added or removed), this should be fast. And the API is really nice now so all this can happen in the background and invisible to the user. I don't think it's a big priority, but if you ever have some runs that you would like to speed up, send me a short example and I'll work on it.

[tigerchenlu98]

Yep very valid about the more general API - certainly if the user is using the C version they can be expected to write out the more verbose functions.

You're right about the sx/sy convention - went ahead and fixed that! I think this was only an issue with TidesSpinEarthMoon.

Let me take a closer look at rebx_spin_potential - yeah it's a holdover from tctl. I assumed it was OK because the Hut & Eggleton frameworks are equivalent, but since the previous implementation didn't include changing spins there's probably something missing - also probably just a good idea to change it to the Eggleton notation for consistency's sake anyway.

Regarding speedups: one idea Greg recently pitched me. Currently, the BS ODE integrator uses the same timestep as the simulation, right? Would it be possible to use different timesteps for the BS and the dynamical integrators, since the timescales associated with the spin precession are so much longer? Since the spin axis typically changes so slowly this shouldn't result in large errors - and since the bsstep is so much more expensive than a WHFast step, should be a very significant speedup too.

[hannorein]

Regarding timesteps. That's certainly possible. However, it will add a whole lot of extra complexity related to making sure the different timesteps are synchronized when needed. I think a lower hanging fruit would be to try and optimize the BS step in that case. For example, if we know the spin changes slowly, we could fix the BS subtimestep, limit the number of iterations to one, and turn off some of the failsafe convergence checks.

[hannorein]

Just a few more words of caution. This is a tricky problem. We have a dynamical system consisting of two parts, the orbits which need some short timestep in order to be resolved, and the spins which change over much longer timescales. Because they are coupled, the shortest timescale of the two parts needs to be resolved. If we don't resolve it, even just for one part, the system is simply not converged. It might work in most cases, but one can always construct a case where things will break (quantitatively and qualitatively differ). For example, imagine a situation where the long BS spin timestep exactly matches the orbital frequency of a planet, or a resonant or secular timescale in the N-body part. Then you get some artificial timestep resonance. It's hard to predict when this might happen. Even without such a resonance, the N-body part will always see the spin part with a time lag because the last datapoint is several N-body timesteps old. That sort of asymmetry already exists right now when the timesteps are equal, but will get much worse if the timesteps differ. If the BS spin timestep is, say 1-10% of the spin timescale (which is where you would want to be for optimal performance), then this error is ~1-10%. This error will accumulate over time and quickly lead to order unity errors.

[dtamayo]

I definitely agree with Hanno that there are way too many things that can go wrong to move away from this as the default, and no obvious way to try to auto-detect the right timescales.

This could be a useful capability for many effects though. I agree that for a BS timestep that is 1-10% of the spin timescale the fractional error is of the same order, but I think thanks to the switching scheme the error wouldn't accumulate right? Aren't we alternating WHFast and BS steps? We analyzed cases like this in the REBOUNDx paper, even under dissipation.

I need to think about this a bit more, but I think we're actually close to this capability making custom operator chains

https://github.com/dtamayo/reboundx/blob/master/ipython_examples/CustomSplittingIntegrationSchemes.ipynb

Couldn't we just add a BS stepper to steppers.c, and then do something like

WHFast(dt)....WHFast(dt)Tides(500dt)Spins(1000dt)Tides(500dt)WHFast(dt)...WHFast(dt)?

The user would set the value of 1000 and check for convergence themselves. This would not work for IAS15 though, and would probably be much more complicated for MERCURIUS.

But let's finish merging the pull request and we can then think about this! So excited for all of these changes. Thank you both so much for all of this work. I am not an expert on tides, but I think for tides_spin_potential we just need to combine Eqs 31 and 37 of the Eggleton paper.

[hannorein]

Yes! A custom splitting scheme is the solution!

[tigerchenlu98]

I agree for sure that it's far too unsafe to use as default! But perhaps a safe mode version for default that can be turned off like WHFast?

@dtamayo just want to make sure I get this right: I can't find the explicit expression for the spin potential defined in calculate_spin_potential anywhere in either the Hut 1981 or Baronett papers - do you know what exactly it's calculating, because it's not entirely clear to me.

The comment for that function is "the conservative portion of tidal potential", which would seem to imply simply equation 31 in Eggleton et. al (the potential of the quadrupole arising from spin & tidal deformation). Combining eq 31 & 37 would give us total energy, including point-particle potential and kinetic energy, which I don't think we want right? Correct me if I'm wrong.

[dtamayo]

You're right, we shouldn't include the first two terms of Eq 37, which are calculated by sim.calculate_energy(). It would be nice to include the last term, so that we can keep track of energy conservation, but you're right again that it's not part of the potential. I think we need to update rebx_spin_potential to just calculate delta Phi in Eq 31, and write a rebx_spin_kinetic_energy to calculate 1/2 I Omega^2.

Not sure if you saw Hanno's pull request thread, but what would you think of me rewriting sx, sy and sz as a reb_vec3d? It seems to me from writing examples like that would be a more convenient interface (which is always for me to tell before writing out a few examples!). If you'd be OK with that, what would you think of using 'Omega' as the spin frequency vector and 'I' as the moment of inertia? I can see why you chose the other names to avoid ambiguity, but I think they might be standard enough that it's OK to reserve them for that?

S is just the spin frequency right? Not the spin angular momentum? There are at least some papers that use S for the spin angular momentum I * Omega (e.g., Fabrycky et al. 2007 on Cassini states with dissipation). Happy to go with what you think best--you've thought more about this!

[tigerchenlu98]

@dtamayo Excellent, just wanted to sanity check that! Both functions are now in - rebx_spin_kinetic_energy isn't actually called anywhere since it's not a part of spin potential. I figure we can either rename rebx_spin_potential to rebx_spin_energy and include both, or have some rebx_calculate_total_energy function if it doesn't already exist. The first option seems more natural to me, but let me know what you think!

Yep, took a peek at @hannorein's pull request - sounds like a great idea, I think it's certainly a more natural way to look at things. So long as it won't impact performance of course.

I think s_x/y/z -> Omega_x/y/z and moi->I are both good ideas, and in addition I'm considering k2->kL, if only to be consistent with the notation in the paper. You're right - for moment of inertia at least I chose moi to avoid confusion with inclination, but it looks like REBOUND uses "inc" mostly so that shouldn't be an issue. If you're planning on restructuring the spin axis feel free to go ahead and pull the trigger on the renamings - I'm happy to do it too.

Yep, s is just spin frequency! We should definitely stick with s = spin frequency just to be consistent with Eggleton's notation in my opinion - but I realize that we could probably be more explicit about that in the paper. Something to change for when we get the review back...

[hannorein]

Performance should be better as there will be fewer param lookups.

[dtamayo]

Awesome! @tigerchenlu98 did I understand correctly that s in the code is Omega in the Eggleton paper (second sentence in your last paragraph is a typo)? I agree, all of this can be changed in the review. I agree that sticking with the notation in Eggleton makes perfect sense. I can change k2->kL. It's true that Omega has more potential for confusion, but I'm pretty sure that is always the standard for spin frequency, not just in the Eggleton paper. I can make those changes.

And I think it makes sense to keep potential and spin_kinetic_energy separate. The idea of having all the potentials separate from a total energy function is so that if you include e.g. gr_potential and tides, you could calculate the right energy by taking the pieces you need. I'll add rebx_spin_kinetic_energy to the python side. Thanks!

[tigerchenlu98]

Yep, s in the code is Omega in the Eggleton paper! Should have written s_x/y/z (now) -> Omega_x/y/z (changed to) is the spin frequency.

[dtamayo]

I just spent an hour going back and forth between papers trying to get straight the definition of the tidal love number. As far as I can tell, what we're calling k2 is Q/(1-Q) in Eggleton Eq 133. Is that right? I didn't find kL or Love number anywhere in the Eggleton paper. You probably worked through this to match things up with the Leconte paper. I think I prefer k2 to kL, partially because that's the notation in Murray and Dermott, and partially because that's the named used in the other tides implementation, so if we use kL it would seem like they're supposed to mean different things. Does that sound OK?

[tigerchenlu98]

Yep, that's right! It is a little confusing because it's never explicitly mentioned in Eggleton et. al and the only direct reference to this relation I could find was in the supplementary text of Naoz 2016. But you can recover that expression by equating the equations of Mardling & Lin with Eggleton's, and I doubled checked it with Rosemary herself.

No complaints about k2, those are valid points too!

[dtamayo]

I've changed sx, sy, sz to a vec_3d Omega (on my newpr91 branch). I'm updating the various examples now. Two quick questions:

1) @hannorein Without a way to do spins[i] = ps[0].params['Omega'] if params returns a vec3d, I've written params to take and return a list for parameters registered as vec3d. Do you see a better way forward?

2) @tigerchenlu98 Last annoying question. Is Eq 133 in the Eggleton paper always valid? If so, would it make for the user to set tau rather than sigma? What value do tides papers usually list? I naively would have expected to set the time lag tau on a constant time lag model. Given that both tau and sigma are just adjustable parameters to capture our ignorance, it seems nice that it's obvious what the units for tau is, and it would mean we would just need one helper function (to go from Q to tau) rather than 2 from Q and tau to sigma.

3) I think we're ready to merge after this. Are we missing anything else?

[hannorein]

1. I think this is a good workaround for now. If anyone answers my numpy bugreport/stackoverflow question, we can update it later.

[tigerchenlu98]

Yep, 133 is always valid, and tau is certainly more typically used. I chose sigma just because that's the most basic value that the equations are actually in terms of - but it would not be difficult to do things in terms of tau instead, and that may be more accessible. I can go ahead and change things around.

Unfortunately updating the syntax seems to have messed with the ml19 migration example... looking into that right now. Hopefully it's a quick fix - I'll let you know when I've figured it out and then we can merge!

[dtamayo]

Hi Tiger, yes! I just fixed one C and one python example to make sure it was working. I can update all the rest of the examples. Just wanted to ask about tau so I could do it all at once. I can do both!

[dtamayo]

I changed from sigma to tau. I think this makes things much nicer for the user! (and computational hit is completely negligible). I left a convenience function to go from Q to tau in C, since it's a pain to calculate the mean motion, but I didn't even add it to the python side. It's longer to write out the function call than to do tau = 1 / (2ps[1].nQ).

I was surprised by the Q values in the Earth and Moon example...I thought the Earth's Q was closer to 10. Is it supposed to be the value at some other time than the present? And should we write out a caveat in the example saying that this would be a case where the approximation is suspect since the Earth's spin isn't synchronized?

I'll finish updating the C examples tonight, and perhaps we can merge tomorrow after you and Hanno have a chance to look at it!

[dtamayo]

@tigerchenlu98 @hannorein I think I've updated all the examples now. I ended up taking out the tau_from_Q helper function in C too...it's really not much longer to write it out in C, and it's probably something the user should realize what they're explicitly doing...I'm a bit confused about it myself. In the paper, we mention that this approximation tau^-1 = 2nQ is only valid for spins synchronized to a circular orbit, but we end up using it in several examples where that is not the case (e.g. spin_kozai). I tried to add notes in the examples--take a look and see if you want to sharpen them.

Unless you see anything else left, I think we're ready to merge! @tigerchenlu98 could you take a look? Thank you both so much for all the work that went into making this project happen, and for your patience on this last piece figuring out a good API. I really think it was worth the effort to get this right! I think a lot of people are going to find this helpful.

[hannorein]

Feel free to merge any time! I'll go over the examples at some point, but not sure when...

[tigerchenlu98]

@dtamayo yes the Q-tau relation was something that gave me quite a few headaches. It's because in the literature equilibrium tide theory is widely used but people only think in terms of Q, since tau has ended up being less intuitive. So in principal while that Q-tau relation is only strictly true for a synchronized circular orbit, in practice people just use it everywhere regardless of the validity. The general consensus seems to be that there's so much uncertainty regarding tidal parameters of exoplanets that the invalidity of this approximation doesn't really matter - but as our knowledge of the tidal dissipation of exoplanets grows hopefully the field as a whole will move towards a more nuanced understanding of Q. For now, the relation is useful for reproducing others' results, and for consistency with Mardling & Lin 2002.

I'll take a look at the Earth/Moon system Q values - I do think you're right in that Q should be closer to 10 for the Earth. It shouldn't change anything qualitative though.

And yes, I think we're almost ready! I just want to make sure the ml19 migration example is up to scratch (to be clear, it's compiling and running fine, so the syntax is all good. The dynamics have just gone weird so I must have just messed something up in the conversions). I'll hopefully figure that out today and I'll let you know when I do - then I think we're good to merge!

[tigerchenlu98]

Hmm it looks like whatever I inadvertently messed up was fixed by your changes @dtamayo - it's outputting exactly as expected now, so crisis averted. If you'd like we can wait the full 2 days to make sure everything's exactly as expected, but I'm happy with how it's running and don't expect anything to differ.

There seems to be a great deal of uncertainty regarding the Earth's tidal Q... Lainey (2016) which we use does seem to be the most up-to-date measurement though. I went through and checked with Q=10 and nothing significant changes, so I'm pretty comfortable keeping the Q we currently have.

The Kozai example was going unstable when tides were added in - don't know how that slipped by me. I just changed it to the paper example, and everything works as expected now. I see you changed the GR implementation to gr_potential - I figured for two structureful bodies the planet would also be worth considering for GR, hence the need for gr_full. I could well be wrong, so happy just use gr_potential if you think that's a good enough approximation.

Your comments look good for the tau-Q relation stuff! It's basically a good order of magnitude starting point for a tidal dissipation parameter, which given the uncertainty around exoplanet dissipation is justified.

I'm good with a merge now!