search

gandalfcode / gandalf

GANDALF (Graphical Astrophysics code for N-body Dynamics And Lagrangian Fluids)
GNU General Public License v2.0
46 stars 12 forks source link

Relative opening #118

Closed rbooth200 closed 7 years ago

rbooth200 commented 7 years ago

I've added a GADGET-2 like relative opening criterion, in the hope of having a force calculation that is accurate even when the terms all nearly cancel.

We need to see whether: 1) This helps 2) The eigenmac criterion isn't already better.

Currently I can't get anything sensible with the Jeans tests at all though...

rbooth200 commented 7 years ago

I believe this should work properly now, althrough there are likely still issues with the Jeans test. I'm still not certain whether PostInitialConditionsSetup is correct, but that's another story.

rbooth200 commented 7 years ago

I've found the bug that was causing the jeans test problem which is now fixed. The relative opening criterion is working well - but the code currently is in need of cleaning up. I'm delyaing it because I want your oppinions on how to do this.

I'm attaching a comparison of the different tree opening methods for the jeans test below. You will see that both the eigen mac and the relative opening criterion do fine, but the geometric version is crap.

Fully opening the tree: i.e. bruteforce gravity Bruteforce

Geometric Geometric

Eigenmac Eigenmac

Relative Relative

giovanni-rosotti commented 7 years ago

Good news! What exactly do you want opnion on? Just in general, or are you thinking to something specific?

giovanni-rosotti commented 7 years ago

Also, can you add a nose test for this?

rbooth200 commented 7 years ago

I'll add a few nose tests. By the way we need to sort this today, since self-gravitating simulations are currently broken in the master branch.

My main questions are: 1) Should we include the relative opening criterion as an option? I guess it fits with the kicthen sink idea, or is the eigenmac sufficient. It's not clear to me that the eigenmac will always work since it's based upon the potential rather than the acceleration, but it does appear to work well for this test case at least.

2) Are you happy with the implementation.

giovanni-rosotti commented 7 years ago

I agree that this is a matter of urgency. I'd say also that we need to make clear in the userguide that the geometric criterion is crap for these problems.

Can you think of problems where the eigenmac criterion fails, and the relative one is good? This would be a strong argument to keep it. I am not familiar with the literature on these things, but I would imagine these are "solved" problems.

I'll now give a look at the implementation.

rbooth200 commented 7 years ago

I can't find a direct comparison between the two. @dhubber mentions them both in the Seren paper, but never tests the Gadget2 one.

giovanni-rosotti commented 7 years ago

Then this is an argument for keeping it, and doing a comparison in the GANDALF paper...

rbooth200 commented 7 years ago

I'm leaning towards keeping it. It would also be more familar for people migrating from Gadget2, which can only be a good thing.

rbooth200 commented 7 years ago

At the moment, we are using two different values for the mac error and the relative opening criterion, along with two string parameters to choose which method to use. This seems like overkill, I was thinking of using a single parameter for each.

Seem ok?

giovanni-rosotti commented 7 years ago

go for it!

rbooth200 commented 7 years ago

How expensive are the multipole moments to compute? They are needed for the eigenmac opening criterion, so it seems to me that we should always calculate them, even when using monopole forces.

giovanni-rosotti commented 7 years ago

"time it" is the only answer i can give...

On Mon, Dec 19, 2016 at 12:28 PM, Richard Booth notifications@github.com wrote:

How expensive are the multipole moments to compute? They are needed for the eigenmac opening criterion, so it seems to me that we should always calculate them, even when using monopole forces.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/gandalfcode/gandalf/pull/118#issuecomment-267952753, or mute the thread https://github.com/notifications/unsubscribe-auth/ABdSgragaalJcRK3FANzNj5RxqU6qbu2ks5rJnhQgaJpZM4LPcpn .

dhubber commented 7 years ago

ok, sorry for being late to this discussion (you started the flurry of messages just as lunch was beginning here).

First of all, it's definitely good to see that this is a very specific problem (i.e. with the geometric opening criterion) and not with other general important things (e.g. tree walk, Ewald). However, I guess the obvious question is, are we sure this is not a bug/programming error instead of a genuine algorithmic weakness? I can see the potential for problems when the net force is quite small/zero but won't this also be true for the other criteria? I might have a peek at the geometric code now just to be sure there's not a bug.

I remember in SEREN I did the relative force criterion using only the gravitational acceleration but I think I removed that so now we only have the total acceleration right? If that's true, then doesn't the tree open more and cells and tend towards brute force the closer to hydrostatic equilibrium is becomes? Maybe this is not a bad thing to have at all since we have higher accuracy when the net force is small (so we are not dominated by noise) but I just want to be sure we understand what is happening. (However, I just noticed there's an 'atree' in the code. Is this what was agrav or is it simply there for debugging purposes?)

Should we include the relative opening criterion as an option? I guess it fits with the kicthen sink idea, or is the eigenmac sufficient. It's not clear to me that the eigenmac will always work since it's based upon the potential rather than the acceleration, but it does appear to work well for this test case at least.

Can you think of problems where the eigenmac criterion fails, and the relative one is good? This would be a strong argument to keep it. I am not familiar with the literature on these things, but I would imagine these are "solved" problems.

Yes, we definitely should include it as an option. It was kind of one the to-do list anyway but never got around to it. The eignemac approach definitely works well in some circumstances and I think it is slightly more conservative than the relative criterion and in particular captures cases where the gravitational acceleration would by chance cancel out but the potential never does so it maintains a better opening criterion. However, there are probably cases where the relative is better but I can't be certain. This would need a broader study of how the trees perform in different density fields (something I guess we don't have time for right now).

How expensive are the multipole moments to compute? They are needed for the eigenmac opening criterion, so it seems to me that we should always calculate them, even when using monopole forces.

They are not super-expensive compared to everything else in the tree build/stock. While it's true they inflate the CPU time for monopole-only, we should really be using quadrupole moments instead of monopoles since the cpu time for a given error is usually lower. We currently only have the faster (Taylor expansion) calculation for monopole and not for the quadrupole, but that's another thing.

rbooth200 commented 7 years ago

Hi David,

Thanks for taking the time to go through this. So I think we are agreed that this a good thing. I believe its (more or less) ready to go. I don't think there is any need to check the tree walk for the geometric case. I was playing around with it and it seems to work as expected. For example, you get brute force as theta -> 0.

I think we'd have to test to see which is generally better. I'd expect that the relative criterion is more robust because it will expand the tree more for those particles with very low gravitational forces. It's worth noting that the gravitational force is used in the tree opening criterion, so it won't go to bruteforce in hydrostatic equilibrium, except if particles at the centre have vanishing acclerations. Probably the eigenmac is a reasonable default as a compromise.

By the way I named the variable atree since it was just the tree acceleration initially, but now I've ended up putting the stars and periodic contribution in. Afterall, if the stars dominate the forces then there is perhaps no need to expand the tree extremely far. But you can come up with counter examples, so I don't think this is ever going to be an exact science... I can rename atree -> agrav if it's highly desired. I could also try to add back the grav / hydro / total acclerations to the diagnostics, but I don't think that's a priority: we should just open an issue.

We currently only have the faster (Taylor expansion) calculation for monopole and not for the quadrupole, but that's another thing.

This is no longer true, I implemented a fast quadrupole sometime in the past...

rbooth200 commented 7 years ago

The current implementaion of the eigenmac is acutally broken since it needs the old value of gpot, which is no longer available when the force calculation is done. Unfortunately this is unavoidable with MPI, without adding another variable to the particle or the tree.

I'm currently getting around this issue for gadget criterion by having the tree store the old values of the acceleration stored in the tree. Any objections to me storing eigenmac factor in the tree instead of amin when using the eigenmac?

I'm no closer to fixing MPI by the way...

dhubber commented 7 years ago

I'm currently getting around this issue for gadget criterion by having the tree store the old values of the acceleration stored in the tree. Any objections to me storing eigenmac factor in the tree instead of amin when using the eigenmac?

No, not at all. Think that's maybe how I did things in SEREN originally anyway!

rbooth200 commented 7 years ago

The zero'ng of gpot is likely causing problems for sink creation too. I'm going to create an issue

rbooth200 commented 7 years ago

Apart from clang timing out, I believe this is now working and can be merged once the tests have finished.

giovanni-rosotti commented 7 years ago

There's only one thing I don't understand: if you have fixed the zeroing of gpot, why is the issue with sinks still a problem?

rbooth200 commented 7 years ago

gpot is stored in the tree cells during the tree build, which is all that is needed for the tree walk. As far as the tree walk is concerned we can then zero everthing before UpdateAllProperties. But doing this means that gpot is not available for flagging potmin, which is done during UpdateAllProperties and is needed for sink-creation.

Solving this problem is non-trivial if iterate the force calculation, since we might zero the potential when we need it later. I'm also not convinced that we are doing the sink creation in an optimal way either. Can this not be done with a tree-walk rather than storing extra variables in the particle and looping over them all? We might need to store gpotmin and rho_max in the tree, but I don't see a problem there.

rbooth200 commented 7 years ago

I think the commit that I've just pushed should solve the problem. When iterating we might compare a new potential with an old one, which means that it's now possible that both a particle and its neighbour will be flagged as potmin, but that was also the case before. So there are at least no new bugs...

But we still need to have a good look at the sinks.

giovanni-rosotti commented 7 years ago

Ok, I think this is enough for this branch then. I'll reference this in the discussion on sinks. As soon as Travis is finished I'll merge this...