Change to Grid_DepositBaryons.C that fixes a vexing problem with the gravitational potential for l=2 grids that are not aligned with the root grid.

#120 Merged
Repository
clean
Branch
week-of-code
Repository
enzo-dev
Branch
week-of-code
Author
  1. Michael Kuhlen
Reviewers
Description

I believe I've identified (and fixed) a bug that was causing a gravitational potential field problem with zoom-in ICs for which the level=2 initial grid is not aligned with the root grid. (This was the root of Tom's PR from a few days ago.)

The following initial grid setup illustrates the problem:

                      TopGridDimensions = 64 64 64
CosmologySimulationNumberOfInitialGrids = 3

    CosmologySimulationGridDimension[1] = 24 26 26
     CosmologySimulationGridLeftEdge[1] = 0.390625 0.390625 0.406250
    CosmologySimulationGridRightEdge[1] = 0.578125 0.593750 0.609375
        CosmologySimulationGridLevel[1] = 1

    CosmologySimulationGridDimension[2] = 32 32 32
     CosmologySimulationGridLeftEdge[2] = 0.421875 0.4296875 0.4453125
    CosmologySimulationGridRightEdge[2] = 0.546875 0.5546875 0.5703125
        CosmologySimulationGridLevel[2] = 2

[MUSIC .conf file (http://pastebin.com/iXGcmwCy) and full Enzo parameter file (http://pastebin.com/9tVh0GPk), if you want to try it out for yourself.]

With the current enzo-dev tip, the potential looks like this after only 1 timestep: http://imgur.com/SpTrg. This weird gradient in the region surrounding the hires region is completely unphysical and causes problematic accelerations.

The problem arises from an incorrect deposition of the l=2 baryons into the l=0 root grid GravitatingMassField (GMF). This is what the l=0 GMF looks like: http://imgur.com/tMiA1. Note the unphysical edge features. [For this plot I turned off particle deposition and the comoving source term, and simply dumped the 70x70x70 GMF to disk immediately after the call to PrepareGravitatingMassField2a() in PrepareDensityField(). The plot is simply an x-slice through the center of the cube.]

The bug occurs in Grid_DepositBaryons.C because dep_grid_cic() expects GridStart >= 0, but when the l=2 grid is offset by half a root grid cell, then GridStart can become negative (-0.5) and this causes dep_grid_cic() to produce bogus results. The problem is that nint(x) (defined in macros_and_parameters.h), which is used to calculate GridOffset, rounds x to the nearest integer, while we actually need greatest-integer-less-than x in this case.

I fixed this by switching from nint() to int() in the calculation of GridOffset, and with this fix the GMF looks as expected (http://imgur.com/pkfzo) and the gravitational potential problem completely disappears (http://imgur.com/752EM).

It seems to me that we may also want to switch from nint() to int()+1 for GridOffsetEnd (imagine a level=3 grid offset by 0.25 root grid cells), but I didn't put that change in since it hasn't caused a problem for me so far. Opinions?

I have not run the test suite with this change, but will hopefully do so over the next few days. If one of you wants to give this a shot, be my guest.

Mike

Update: I incorporated a change suggested by Greg Bryan, see the comment below.

Update 2: I had forgotten to remove the -1 from the updated GridOffset calculation, so effectively I has replaced nint(x)-1 with int(x)-1, which doesn't make any sense. This update fixes it to the intended int(x).

Update 3: Siwtched back to nint(x)-1. The only change in this PR is now the switch from CIC to NGP for subgrid baryon deposition. This solves the subgrid potential problem that was the root of this PR.

Comments (14)

  1. Michael Kuhlen author

    I finally got around to running the quick test suite on this change. I started from the enzo-2.2 changeset (301413545861), which is the basis for the current Enzo gold standard on Amazon.

    With the changes in this PR (nint -> int), all tests pass except for GravityTestSphere (21 failures) and PhotonShadowing (1 failure): test_results.txt.

    I looked at GravityTestSphere using Sam's scripts that he wrote for the GravityTestSphere discussion in PR #114. As you can see from this plot, this PR's change goes in the same direction as the effect of Greg's revert from PR #115, but not quite as far. Including both changes gives identical results to Greg's fix.

    I haven't investigated the PhotonShadowing failure, but since that is one of the tests that also fails for the current tip (due to code changes since the 2.2 release), I'm not going to worry about it.

    In summary, I believe that once the gold standard is updated for the most recent tip, the changes in this PR will pass all quick tests.

  2. Greg Bryan

    Hi Mike -- I really appreciate you digging into this and finding the problem! I'm still trying to understand exactly what is going on, so can we hold off on accepting for a little while? Your results look great, so I'm very tempted just to hit pull, but I'm missing something, probably very simple.

    My issue is that I don't quite understand your explanation of what is going wrong -- it makes sense, but if I say x = (GridLeftEdge[dim]-TargetGrid->GravitatingMassFieldLeftEdge[dim])/TargetGrid->GravitatingMassFieldCellSize (just for notational ease), then (in the current code, after my change, but before yours):

    GridOffset = nint(x) - 1

    and

    GridStart = x - GridOffset

    or GridStart = x - nint(x) + 1

    in which case I don't see how GridStart can become less than 0. I think that x - nint(x) should range from -0.5 to 0.5, and so GridStart should range from 0.5 to 1.5, and should never be negative. What am I missing?

    Thanks, Greg

  3. Michael Kuhlen author

    I've looked into this some more, and Greg and I have exchanged a few emails about the issue. I agree with Greg that with his fix in PR #115, it's not possible to get a negative GridStart and hence dep_grid_cic() isn't ever passed a negative value. So that was apparently a bit of red herring that only applied to the pre-PR#115 state of Enzo.

    However, it turns out that even with Greg's fix the potential in the cosmological example given above still isn't quite right, see here.

    As explained in more detail in an (upcoming) email on enzo-dev, this turns out to be due to using the CIC deposit scheme for the refined grids, which in this example introduces a slight density fluctuation along the subgrid edges. In normal use cases this effect would be smaller than the cosmological initial density perturbations and hence wouldn't be noticed. But for very small initial perturbations, like in this example, it becomes a problem. The solution may be to switch to nearest neighbor interpolation for the subgrid deposition only, and Greg has implemented this change. It gives very good results, both for this cosmological test problem and for TestGravitySphere.

    I've ported these changes to this PR, which is now a combination of switching from nint()-1 to int() and switching to NGP for subgrid deposition. This is not yet the final word, and it may indeed be better to keep nint()-1 even with NGP, so please don't accept yet.

  4. Greg Bryan

    Mike, thanks for summarizing and posting the exchange on enzo-dev. Thanks also for updating this PR!

    My only comment is that I think it would be better to keep the nint(x) -1 for the first computation of GridOffset (lines 124-128 in the modified version) if doing CIC interpolation (and use int(x) if doing NGP interpolation). The nint(x) -1 really is correct for CIC.

    Also, I think the second computation for GridOffset (for the case of TargetGrid == this, on line 131 of the modified version) should be nint. This is because in this case we don't want to go beyond the active region and the computed number should really be an integer (since the difference computed is an integer number of CellSizes) -- the nint in this case just makes sure we don't convert 3.99999 to 3.

  5. Michael Kuhlen author

    Hi Greg, thanks for your comment. I have a few more questions. Some of these may be a bit naive and simply reflect a lack of familiarity with the inner workings of this code.

    1) I'm unclear about why we should still use CIC for the root grid. If we're depositing baryons from grids on the same level, then why is there any need for interpolation at all? Shouldn't the mass in a given l=0 cell be just that, and not spread out over multiple cells of the target grid? The code would become easier to understand if we could just use NGP all the time...

    2) Is the only reason for using nint() to avoid the 3.99999 -> 3 situation? Because of course this situation could also occur for subgrids, when they happen to be aligned with the root grid... So maybe instead we could always use int() (even when TargetGrid==this), but then check afterwards whether this undesired truncation has occurred. So I guess something like this:

    float GridOffsetFloat = (GridLeftEdge[dim] - TargetGrid->GravitatingMassFieldLeftEdge[dim]) /
                            TargetGrid->GravitatingMassFieldCellSize;
    GridOffset = int(GridOffsetFloat);
    if ( (GridOffset + 1.0 - GridOffsetFloat) < 1e-6 )   // or whatever tolerance is appropriate
      GridOffset += 1;
    

    I think this should always give the first root grid cell (lowest index) to contain any part of the grid to be deposited. GridOffset may need to lowered by 1 if we're using CIC.

    3) I'm a little confused about this part:

    if (TargetGrid == this)
        GridOffset[dim] = max(GridOffset[dim],
                              nint((TargetGrid->GridLeftEdge[dim] - TargetGrid->GravitatingMassFieldLeftEdge[dim]) /
                                   TargetGrid->GravitatingMassFieldCellSize) );
    

    If TargetGrid == this, then TargetGrid->GridLeftEdge is identical to GridLeftEdge, and in that case the second argument of the max(,) call is just GridOffset[dim] + 1, and so it will always return exactly that. So we might as well write, 

    if (TargetGrid == this) GridOffset[dim] += 1;
    

    No?

  6. Greg Bryan

    Hi Mike -- These are excellent questions and I'm not sure I have excellent answers, but here is my current thinking.

    1) I'm unclear about why we should still use CIC for the root grid. If we're depositing baryons from grids on the same level, then why is there any need for interpolation at all? Shouldn't the mass in a given l=0 cell be just that, and not spread out over multiple cells of the target grid? The code would become easier to understand if we could just use NGP all the time...

    Well, The point of using CIC deposition is two-fold: (1) in principle lower noise, although for the top-grid this is not really an issue, but more importantly (2) to allow us to make an estimate of the time-centered baryon field. The grid positions are shifted slightly (by velocity * dt/2) before depositing on to the grid. This should provide us with an easy estimate of the rho(n+1/2), which is then used to compute the accelerations (at time n+1/2), which is assumed by the code in UpdateParticlePositions. Now this is probably not the best way to time-center the gravitational acceleration -- in fact lately, I've been playing around with some ideas for changing this (one way would be to make an estimate of the advection by a 1/2 step)., so maybe we should switch entirely to CIC; however, that would be a larger algorithmic change than this PR is currently about, and would probably require more testing. Note that none of the tests we've done so far have tested the time-dependence!

    2) Is the only reason for using nint() to avoid the 3.99999 -> 3 situation? Because of course this situation could also occur for subgrids, when they happen to be aligned with the root grid... So maybe instead we could always use int() (even when TargetGrid==this), but then check afterwards whether this undesired truncation has occurred. So I guess something like this:

       float GridOffsetFloat = (GridLeftEdge[dim] - TargetGrid->GravitatingMassFieldLeftEdge[dim]) /
                             TargetGrid->GravitatingMassFieldCellSize;
       GridOffset = int(GridOffsetFloat);
        if ( (GridOffset + 1.0 - GridOffsetFloat) < 1e-6 )   // or whatever tolerance is appropriate
           GridOffset += 1;
    

    I think this should always give the first root grid cell (lowest index) to contain any part of the grid to be deposited. GridOffset may need to lowered by 1 if we're using CIC.

    For the computation of GridOffset when using CIC:

        GridOffset[dim] = nint((GridLeftEdge[dim] -
                    TargetGrid->GravitatingMassFieldLeftEdge[dim])/
                   TargetGrid->GravitatingMassFieldCellSize) - 1; 
    

    the nint is used for a different purpose -- here we are making use of the 0.5 in nint(x) = int(x + 0.5) because we know that the CIC deposition will extend 0.5 cell widths from the position of the particle. Maybe this line would be clearer if we actually used int(x + 0.5) instead of nint to make it conceptually clearer that there is a factor of 0.5 in this calculation? [In fact, it would probably be even clearer if we used int(x - 0.5) instead of int(x + 0.5) - 1 ].

    For the NGP calculation, we don't need the 0.5 part because the NGP deposition width is 0.0, so here we use int(x).

    We could do the tolerance check you suggest -- for the NGP calculation this would ensure that we never use a region which is larger than necessary, although note that if we do, there is no error, we just pass around slightly more data than necessary. I have shied away from this because of the difficulty in setting the tolerance -- for example, in this case I would use 0.5*CellWidth[dim]/TargetGrid->CellWidth[dim] instead of 1e-6.

    For the CIC, your suggestion also would work (making the tolerance check and then subtracting 1), although it does prescribe a larger region than necessary if this != TargetGrid (although our suggested change makes this case always into NGP so perhaps that's moot, but I would prefer not to make assumptions like that).

    3) I'm a little confused about this part:

       if (TargetGrid == this)
          GridOffset[dim] = max(GridOffset[dim],
                                 nint((TargetGrid->GridLeftEdge[dim] - TargetGrid->GravitatingMassFieldLeftEdge[dim]) /
                                    TargetGrid->GravitatingMassFieldCellSize) );
    

    If TargetGrid == this, then TargetGrid->GridLeftEdge is identical to GridLeftEdge, and in that case the second argument of the max(,) call is just GridOffset[dim] + 1, and so it will always return exactly that. So we might as well write, 

        if (TargetGrid == this) GridOffset[dim] += 1;
    

    No?

    Well, yes, that's true, but I wanted this statement not to depend on what was done with GridOffset in the first part of the calculation. As we have seen, there is some uncertainty as to exactly what the correct value of GridOffset should be (and when I was developing this code, I messed around with different values, a bit like we're doing now). The point of this statement was to make sure that we don't try to deposit outside of the ActiveRegion of this grid, regardless of how GridOffset was computed, since that would be an actual error (rather than just passing around extra zeros).

  7. Michael Kuhlen author

    Thanks for those explanations, Greg. Given the fact that using a region that is larger than necessary doesn't give the wrong answer, I think it's probably simplest and cleanest to just stick with the original nint(x)-1 in all cases, even for NGP. This minimizes the code changes, and it gives the desired results for this cosmological test problem as well as for TestGravitySphere. I have updated the PR, and I think this PR is ready to be accepted. I believe we can also close (decline) PR#114.

  8. Nathan Goldbaum

    Hi all,

    This shouldn't hold up the acceptance of the PR but I'm curious if we understand why enzo misses the analytic solution for the potential in the TestGravitySphere problem. I understand that the force error is very small but won't there be issues with energy conservation if, say, a particle falls to the bottom of the potential and picks up a bit of extra spurious energy due to the error in the potential?

  9. Greg Bryan

    I think this PR should be merged. I've been using it for a while and it seems stable. It also significantly improves the solution to the Bertschinger self-similar test (see discussion in PR#161)

    In response to Nathan's question above, the problem is unrelated to this PR and was fixed in PR#124.

  10. Nathan Goldbaum

    This PR is necessary to reproduce the analytic result for the SelfSimilarInfall problem. Not only that, but looking over the code and the PR comments, it seems that this is a sensible, well thought out change that will improve answers for zoom-ins and hierarchies that have offsets between the level 2 grid and the root grid.