In this fix I propose only correct baryons to the mid time step for the current level and one level of child grids rather than the current way of trying to deposit all baryons from all finer levels. This now avoids errors in mass conservation that could occur when grandchildren grids were not aligned with their grandparent. Sam Skillman has helped a ton in testing this already on AMRCosmology and TestGravitySphere next to the NestedCosmologRun I used to diagnose it. Acknowledgements are also due to John Wise and Mike Kuhlen who noticed the problem in the first place and also narrowed down where the bug may be hiding in the first place.
We discovered that Tom's changes lead to the correct, analytic solution in the region outside the collapsing cloud. This is an improvement over the current situation, which misses badly in the region outside the cloud.
Thanks Nathan for summarizing! Let me add another 2 cents
The potential may not be the most important thing to look at …
When I differentiate the radial potential we measured to get the radial force my conclusion is quite different then just from looking at the potential.
I compare the old, the reference(old with much more resolution and buffer zones), the analytic and the new solution in the force.
There is little change in the forces near the cloud but very large differences at large scales. The force plot very clear would argue for the new DepositBaryons.C
I include the force plot and the gravitational potential one below.
I think our next tests should have a dynamical collapse. Given the force plots I'd expect the new version to work just as well as the old one for the central parts of spheres but be a significant improvement on the outside.
May the force be with us,
PS: thanks all the folks on irc today. It was a great enzo day with everyone chipping in, people working together, discussing, learning together. Very fun!
There is indeed a bug with the current code, but I'm not sure this is the correct solution to the problem. Instead it would be better to work out why the deposition is not working in the first place.
I think The key observation is that the potential on the root grid is completely flat -- that's telling us that the gravitating mass field on the root grid is completely featureless, which in turn means that no mass from the subgrids is getting deposited there. After a little digging, I found that, indeed, this is the case -- Grid_DepositBaryons is not calling the CIC deposit routine on any grids more than the 2 levels finer. The problem turns out to be due to some changes from about a year ago (I think), that changed how Grid_DepositBaryons.C computes the GridOffset and GridOffsetEnd of the little temporary grid that holds the deposited baryons. Because of the error in the computation of the size, this covering grid is sometimes too small. There always used to be a check for this which generated an error; however, at the same time as these changes were made, this check was deactivated so that the routine just returns if the covering grid is too small.
If I revert these changes, I get a potential result which is very similar to Tom's changes (http://www.astro.columbia.edu/~gbryan/gpot/gpot.png), and so presumably has the same good match in the acceleration that Tom found. (By the way, I suspect the reason the acceleration matches, while the potential does not is a slight potential mis-match at boundaries which is not reflected in the acceleration).
There is an additional question, which is: is it necessary to deposit baryon subgrids, or can we simply deposit a few levels. This is an approximation, but may be a good one in many cases. I think it could be useful to add that, but it should be done as a runtime parameter, and I think that it should default to the current behaviour (so somewhat different from the way this is coded here).
When I was trying to debug a deeply nested grid simulation, I reverted Tom's +/-1 change to the grid boundaries in the deposit routine. With this reversion, I still saw the artifacts (i.e. missing mass) in the potential field when the nested grids weren't aligned with their grandparent grids. I still think there's some bug within the deposit routine that's causing these artifacts because, in principle, it should work when depositing any descendant grid to the coarse grid, and I tend to agree with Greg on this one.
Ah indeed. It surely is great that we have pull requests now. What I was trying to do over a year ago was to use gradients in the gradients of the potential for refinement. div grad Phi would just give the density. However, the mixed derivatives allow one also define something very analogous to a Jeans criterion that would ensure proper resolution for the baryons even when gravity is dominated by the tidal field of some nearby heavy object. This occurs quite frequently, e.g. in accretion disks in star formation, black holes, and also in galaxy formation simulations when stars and DM eventually dominate the potential and it is not the Jeans criterion counting only the self-gravity of the gas which is relevant. I kept having the strangest problems with it though and kept trying to sort this out. This lead to that awful commit which introduced the bug! I had tested it on 3 different problems but all of them were symmetric around the midpoint and did not show the problem we are finding now ... It is great we have pull requests and a testing framework that will allow all of us to make sure I don't introduce similarly bad bugs again!
I totally agree with Greg, that a runtime parameter LevelsUpToDepositAdvectedBaryons would make the most sense. I suspect that the problem that I could not get the TotalJeansLength refinement to work was because the extrapolated "positions" of the baryons led to an implausible density. With this runtime parameter it could help in keeping the predicted density at the mid time-step to be smooth on the scale of the grid under consideration.
Ok. We keep having an issue though.
Reverting as Greg suggested definitely works for fine for the sphere test and give us the same results. However, it still has somewhat of a trouble in a Cosmology case. The change I introduced a year ago allows potentials to be contiguous in cosmology case:
while the old way gives:
These offsets would normally indicate that there is some difference in masses (surface integral of forces around grid patch = two integral of density over the patch volume).
However, the change I made long time ago then broke what happens with grand-children and my symmetric tests all missed that. The "fix" of this pull request then tried to sort this out again and does give in our current tests identical answers for the TestGravitySphere and also still looks better for the cosmology run.
So we will need to do more tests to see what the need is to do DepositBaryons recursively all the way down to the finest level. If we can find a test case that shows this is in fact important than perhaps we revert to the original enzo way otherwise there is a chance we could just use what we have in this pull request.
Does anyone have a suggestion of what such a test may be? It likely would have to be a dynamic test. Perhaps an advecting sphere in hydrostatic equilibrium to see something change.
Perhaps we may conclude that the full deposit is at least in these not to asymmetric cases not necessary.
I better try this one now displacing the sphere slightly from the center and look for more potentially interesting cases that improve our test-suite to define the metric of what any change to the GravitySolver has to pass.
Anymore ideas along these lines would be appreciated.