Symetrical FMM solver that conserve linear momentum#699
Symetrical FMM solver that conserve linear momentum#699Yrisch merged 30 commits intodanieljprice:masterfrom
Conversation
…on tested on binary setup)
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
There was a problem hiding this comment.
Summary of Changes
Hello @Yrisch, I'm Gemini Code Assist1! I'm currently reviewing this pull request and will post my feedback shortly. In the meantime, here's a summary to help you and other reviewers quickly get up to speed!
I've reworked the Phantom self-gravity solver to introduce a symmetrical Fast Multipole Method (FMM). The primary goal of this update is to ensure robust conservation of linear momentum, which is crucial for the accuracy and stability of N-body simulations. This new approach employs a dual tree walk and a refined Taylor expansion to calculate gravitational interactions. While maintaining performance comparable to the previous method, this implementation provides a significant improvement in momentum conservation, validated by new dedicated tests.
Highlights
- Symmetrical FMM Solver: I have implemented a new symmetrical Fast Multipole Method (FMM) solver for self-gravity, which significantly improves the conservation of linear momentum in simulations. This addresses a key challenge in N-body problems by ensuring the total force on the system remains zero.
- Dual Tree Walk Implementation: The new solver utilizes a dual tree walk algorithm, similar to approaches by Dehnen (2000) and Springel (2021). This involves traversing two k-d trees simultaneously to efficiently compute node-node interactions and propagate forces through Taylor expansions.
- Revised Multipole Moment Formulation: I have updated the formulation of Cartesian multipole moments (order=3) and Green function derivatives. This includes changes to how quadrupole moments are calculated and how force and potential contributions are propagated and combined, ensuring accuracy and momentum conservation.
- Enhanced Momentum Conservation Testing: New unit tests have been added to
test_gravity.f90specifically to verify the linear momentum conservation of the FMM solver. These tests confirm that the sum of all forces in the system is conserved to machine precision, particularly in binary setups.
Using Gemini Code Assist
The full guide for Gemini Code Assist can be found on our documentation page, here are some quick tips.
Invoking Gemini
You can request assistance from Gemini at any point in your pull request via creating an issue comment (i.e. comment on the pull request page) using either /gemini <command> or @gemini-code-assist <command>. Below is a summary of the supported commands.
| Feature | Command | Description |
|---|---|---|
| Code Review | /gemini review |
Performs a code review for the current pull request in its current state. |
| Pull Request Summary | /gemini summary |
Provides a summary of the current pull request in its current state. |
| Comment | @gemini-code-assist | Responds in comments when explicitly tagged, both in issue comments and review comments. |
| Help | /gemini help |
Displays a list of available commands. |
Customization
To customize Gemini Code Assist for GitHub experience, repository maintainers can create a configuration file and/or provide a custom code review style guide (such as PEP-8 for Python) by creating and adding files to a .gemini/ folder in the base of the repository. Detailed instructions can be found here.
Limitations & Feedback
Gemini Code Assist may make mistakes. Please leave feedback on any instances where its feedback is incorrect or counter productive. You can react with 👍 and 👎 on @gemini-code-assist comments or fill out our survey to provide feedback.
You can also get AI-powered code generation, chat, as well as code reviews directly in the IDE at no cost with the Gemini Code Assist IDE Extension.
Footnotes
-
Review the Privacy Notices, Generative AI Prohibited Use Policy, Terms of Service, and learn how to configure Gemini Code Assist in GitHub here. Gemini can make mistakes, so double check it and use code with caution. ↩
There was a problem hiding this comment.
Code Review
This pull request refactors the self-gravity solver to implement a symmetrical Fast Multipole Method (FMM) aimed at better linear momentum conservation. The changes are extensive, introducing a dual tree walk and a new formulation for multipole moments and their propagation via Taylor series.
My review has identified a few issues. The most critical one is in the new compute_M2L routine, where the quadrupole contributions to the force derivatives seem to be missing, which could impact the accuracy of the gravity calculation. This is further compounded by the removal of the finite difference tests that would have verified these derivatives. I've also found a minor issue with an error message and a tautological check in the new tests.
Overall, this is a significant physics improvement, but the potential issue with the force derivative calculation should be addressed to ensure the correctness of the new solver.
|
In the new formulation, the moment contributions are only applied to the force terms. No warnings here. |
|
@themikelau would be great if you could try this branch on a common envelope simulation, should improve linear momentum conservation, but would be good to check that there are no other significant issues in the gravity calculation |
|
seems like a genuine error trying to run SETUP=converging here... |
|
Yes, there are still some bugs on this branch. MPI is broken. My solution fell apart. I thought that every process could compute the long ranges without comm. But obviously it was too good to be true. It will need more sophisticated algorithms. |
|
The converging setup is fixed. I disabled the sym FMM with MPI until we have a solution. The code will use the original algorithm if launched with MPI. |
|
We lose the 40% because of |
|
@themikelau we believe the major slowdown here is now fixed, please can you check [a portion of] the calculation again from your side to clarify the performance hit? We are thinking any remaining performance loss may be resolvable by taking a slightly more generous opening criterion... |
|
I can do that but I am quite busy this week. Do you expect worse linear momentum conservation compared to my last test? |
|
Just to clarify, is this branch up-todate (the last commit was two weeks ago)? |
|
Yes it should ! |
no massive rush, linear momentum conservation should be similar, but the performance should be much improved |
|
Okay. I include updated plots that now have "FMM" (latest changes) and "FMM old" (the previous version that is 5 times slower). My summary is that everything works and the latest implementation is only 1.25 times slower than the control case, while linear momentum conservation is on the level of 10^-13. Comparison of density slices in the orbital plane. Looks alright here. Comparison of x-, y-, and z-components of the centre of mass: Total momentum: Angular momentum: Fractional errors in energy components: Separation between donor's core and companion: Evolution in bound envelope mass: Fractional error in bound envelope mass: |
|
Thanks a lot @themikelau ! 1.25 is much more manageable than 5! 😅 |
|
I would consider paying 25% more for this precision. But getting it down to only ~10% more would be much better :) |
|
Yes, agree! It should be doable... |
|
I think we should just increase the default tree opening a bit to gain back the 25%, then merge! |
|
@Yrisch If this is ready, I can give it another check with a common-envelope simulation. |
|
Yes! Ready to retry. |
|
test failure looks like tolerance just needs increasing slightly (assuming nothing changed here): |
|
Yes nothing changed... |
|
Running the test now... |
|
I have added a comparison with results from the latest branch where the tree opening criterion is increased from 0.5 to 0.55. The performance is much more comparable, with the FMM only being ≈ 5% slower than the control. The increase in tree opening criterion did not significantly affect the level of drift in the centre of mass. Total momentum conservation is slightly worse but still excellent compared to before. I am happy with the branch to be merged. Comparison of density slices in the orbital plane. No significant differences observed. Comparison of x-, y-, and z-components of the centre of mass: Total momentum: Angular momentum: Fractional errors in energy components: Separation between donor's core and companion: Evolution in bound envelope mass: Fractional error in bound envelope mass: |
|
Awesome! Let's merge! |
3 participants
Description:
This is a rework of the Phantom self-gravity solver, implementing a symmetrical fast multipole method that ensures a good conservation of the linear momentum. The method followed is similar to Dehnen 2000 and Springel 2021. It uses a dual tree walk to compute node-node interactions and Taylor expand between nodes to propagate these to leaf nodes.
Like the previous method, it uses Cartesian multipole moments (order=3) and Green function derivatives to estimate the potential and forces in the system. It is worth noting that the formulation of these moments is different from the previous one.
The dual tree walk implementation is really similar to the previous method, giving similar perfs. However, the improvement in momentum conservation is significant up to machine precision. ( in a binary setup)
Components modified:
Type of change:
Testing:
In
test_gravity.F90, we compare the self-gravity solution with a direct algorithm. The results are in tolerances.We also check if the sum of all forces is zero (linear momentum conservation). We test the node interactions with two widely separated blobs of gas and check if the same as the previous tests. We also tested the linear momentum conservation in time by computing a few orbits of a binary (blobs of gas) using
SETUP=binaryDid you run the bots? yes
Did you update relevant documentation in the docs directory? no
Did you add comments such that the purpose of the code is understandable? yes
Is there a unit test that could be added for this feature/bug? no