Correction in the computation of the gravitational potential RHS for the MultiGrid solver

#124 Merged at 99d895b
  1. Greg Bryan

This fixes a long-standing problem in the computation of the RHS of the potential equation in Grid_SolveForPotential.C. This was pointed out by Dan Reynolds last year, but I didn't work out the ramification and the fix until now. The basic problem is that Grid_SolveForPotential.C assumes the cell spacing on the grid is L/Nx, while the multigrid solver assumes it is L/(Nx-1). As Dan points out, this means that we end up solving a problem of the form Laplace(Potential) = (G/a)(GravitatingMass)(dim1-1)(dim2-1)(dim3-1)/dim1/dim2/dim3 as opposed to: Laplace(Potential) = (G/a)*(GravitatingMass)

This ends up solving a persistent error in the potential seen most clearly in the TestGravitySphere test -- the potential is off by typically a value of about 6-8%, entirely due to this error. Using Sam's scripts from [here] ( as can be seen in this image:

I spent some time figuring out to change the multigrid solver, but that turned out to be a bad idea, since it works best with that cell spacing. Instead, it is much easier to adjust the RHS computation. Note these changes are completely separate from the recent discussion about depositing subgrid cells.

This will require a regeneration of the gold standard.

Comments (3)

  1. Sam Skillman

    Hi Greg,

    This seems like an automatic accept to me. Any reason why it shouldn't be accepted, and have the gold standard regenerated?

  2. Tom Abel

    Fantastic sleuthing, Greg! It also gives a noticeable improvement at the grid boundaries. I'll try my TotalJeansRefinement ideas again now. This could have been what was missing. Thanks!