SetAccelerationBoundary (SAB)¶
One of the minor bugs in Enzo that was uncovered by the addition of MHD-CT is the boundary on the gravitational acceleration field.
Enzo currently solves gravity in two phases: first by Fast Fourier Transform on the root grid, then by multigrid relaxation on the subgrids. Unfortunately, each subgrid is solved as an individual problem, and is not very concious of its neighbours.
The problem with this is the ghost zones. Enzo MHD-CT is not a
divergence free method, but a divergence preserving method.
There isn’t a mechanism that reduces the divergence of the magnetic
field. Unfortunately, inconsistencies in any fluid quantity can
lead to divergence in the magnetic field. The magnetic field is
stored on the faces of each computational zone, and are updated by
an electric field that is stored on the edges. Since this data sits
in the face of the zone, whenever two grids abut, they share a
face, so it is vital that both grids describe everything in the
stencil of the face centered fields identically, otherwise they
will get different results for the magnetic field on that face, and
divergence will be generated. It was noticed that in the case of
the AccelerationField
that due to the isolated nature of the
gravity solver, the ghost zones of a subgrid didn’t necessarily
equal the active zones of grids that were next to it. Thus the
Magnetic fields in the shared face would ultimately be computed
slightly differently, and divergence would show up.
The proper fix for this is replacing the gravity solver with one that is aware of the entire subgrid hierarchy at once, but this is quite costly in both programmer time and in compute time. Work has begun on this project at the LCA, but has not yet been finished.
As an intermediate step, Enzo was hacked a little bit. Initially,
the main loop in EvolveLevel.C
looked like this:
for( grid=0, grid< NumberOfGrids, grid++){
Grid[grid]->SolvePotential
Grid[grid]->SolveHydroEquations
}
Among, of course, many other physics and support routines. This was
broken into two loops, and a call to SetBoundaryConditions()
as
inserted between the two.
for( grid=0, grid< NumberOfGrids, grid++){
Grid[grid]->SolvePotential
}
SetBoundaryConditions
for( grid=0, grid< NumberOfGrids, grid++){
Grid[grid]->SolveHydroEquations
}
However, since SetBoundaryConditions()
doesn’t natively know about
the AccelerationField
, another kludge was done. A new set of
pointers ActualBaryonField
was added to Grid.h
, and the true
pointers are saved here, while the BaryonField
array is temporarily
pointed to AccelerationField
. This saved a substantial rewrite of
the boundary setting routines, at the expense of some
less-than-ideal code.
This is not a bug that makes much difference overall in cosmology
simulations, and it does not solve the problem of artificial
fragmentation that has been noticed by some groups. Cosmology tests
have been done that compare solutions both with and without this
fix, and only negligible changes appear. So for most runs, it
simply adds the expense of an extra boundary condition set.
However, with MHD-CT runs it is absolutely necessary, for explosive
divergence will show up. Additionally, and other simulations that
are extremely sensitive to overall conservation or consistency will require
this flag. In any condition where the user is potentially concerned about
we suggest running a test both with and without SAB
, and comparing the answers.
SAB
brings the compuational expense of an additional boundary condition call, and
the memory expense of three global fields, since without it the AccelerationField
exists
only on a single grid at a time, while with it all three fields must be created on the entire hierarchy
at once. This is not a major expense on either count for most simulations.
This is controled by the preprocessor directive SAB
. If this is
defined, the necessary steps are taken to call the acceleration
boundary. In the file machine make file, Make.mach.machine-name
, this should be
added to the variable MACH_DEFINES