The Physics andAlgorithms of Enzo
John Wise (Georgia Tech)!

Most of the slides are originally written by Greg Bryan (Columbia)

Enzo Workshop – 21 November 2014

enzo-project.org

Introducing• Block-structured AMR + N-body!
• Proper or comoving coordinates!
• N-body: Adaptive particle-mesh solver!
• (Magneto-)Hydrodynamics: !
• High-resolution shock capturing scheme!
• OR Finite differencing!
• Chemical network solver (H, He, H , HD)!
• Star & BH formation and feedback!
• Radiative transfer!
• Adaptive (angular) ray tracing!
• OR Flux-limited diffusion2

I. Hydrodynamics

Fluid Equations - grid::SolveHydroEquations
Mass conservation
Momentumconservation
Energy conservation

Ideal Gas EOS
Self-gravity

Field names: Density, Pressure, TotalEnergy, InternalEnergy,
Velocity1, Velocity2, Velocity3

grid class: accessing the fields – grid.h In grid class:BaryonFields[] – array of pointers to each field Fortran (row-major) ordering within each field
GridRank – dimensionality of problem



GridDimensions[] – dimensions of this grid

GridStartIndex[] – Index of first “active” cell (usually 3)
 First (and last) three cells are ghost or boundary zones

int DensNum = FindField(Density, FieldType, NumberOfBaryonFields);
int Vel1Num = FindField(Velocity1, FieldType, NumberOfBaryonFields);
for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) {for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) {for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++) {BaryonField[Vel1Num][GINDEX(i,j,k)] *= BaryonField[DensNum][GINDEX(I,j,k)];}}}Enzo file name convention
 General C++ routines:

Routine name: EvolveLevel(…)
 In file: EvolveLevel.C
 One routine per file!
 grid methods:

Routine name: grid::MyName(…)
 In file: Grid_MyName.C
 Fortran routines:

Routine name: intvar(…)
 In file: intvar.src .src is used because routine is fed first through C preprocessor

PPM Solver: grid::SolvePPM_DE
 HydroMethod= 0 PPM: e.g. mass conservation equation


Flux conservative form:
Mass flux across j+1/2 boundaryIn discrete form:

How to compute mass flux?

 Note: multi-dimensions handled by operating splitting

grid::xEulerSweep.C, grid::yEulerSweep.C,grid::zEulerSweep.C

Grid::SolvePPM_DE// Update in x-directionfor (k = 0; k < GridDimension[2]; k++) {if (this->xEulerSweep(k, NumberOfSubgrids, SubgridFluxes,
GridGlobalStart, CellWidthTemp, GravityOn,
NumberOfColours, colnum) == FAIL) {fprintf(stderr, "Error in xEulerSweep. k = %d\n", k);
ENZO_FAIL("");}} // ENDFOR k// Update in y-directionfor (i = 0; i < GridDimension[0]; i++) {if (this->yEulerSweep(i, NumberOfSubgrids, SubgridFluxes,GridGlobalStart, CellWidthTemp, GravityOn,
NumberOfColours, colnum) == FAIL) {fprintf(stderr, "Error in yEulerSweep. i = %d\n", i);
ENZO_FAIL("");}} // ENDFOR iPPM: 1D hydro update: grid::xEulerSweep
 Copy 2D slice out of cube
 Compute pressure on slice (pgas2d)
 Calculate diffusion/steepening coefficients (calcdiss) Compute Left and Right states on each cell edge (inteuler) Solve Reimann problem at each cell edge (twoshock) Compute fluxes of conserved quantities at each cell edge(euler) Save fluxes for future use
 Return slice to cube

PPM: reconstruction: inteuler
 Piecewise parabolic representation:
qLqqR

 Coefficients (Δq and q6) computed with mean q and qL, qR.
 For smooth flow (like shown above), this is fine, but can cause aproblem for discontinuities (e.g. shocks)
 qL, qR are modified to ensure monotonicity (no new extrema)PLM: reconstruction• Piecewise linear method!
• More diffusive reconstruction scheme, butmore stable.

PPM: Godunov method: twoshock To compute flux at cell boundary, take two initial constantstates and then solve Riemann problem at interface
left stateright state

rarefaction wave
contact discontinuity
shock Given solution, can compute flux across boundary Advantage: correctly satisfies jump conditions for shockPPM: Godunov method: inteuler, twoshock For PPM, compute left and right states by averaging overcharacteristic region (causal region for time step Δt)

 Average left and right regions become constant regions tobe feed into Riemann solver (twoshock).

PPM: Eulerian corrections: euler
 Eulerian case more complicated because cell edge is fixed.
Characteristic region for fixed cell more complicated:

SUBSONIC CASE
SUPERSONIC CASE

Note that solution is not known ahead of time so two-stepprocedure is used (see Collela & Woodward 1984 for details)

Other Riemann solvers• HLL: (Harten-Lax-Leer)!
!
!
!

• HLLC: HLL but considering the contactwave

Difficulty with very high Mach flows
 PPM is flux conservative so natural variables are mass,momentum, total energy
 Internal energy (e) computed from total energy (E): Problem can arise in very high Mach flows when E >> e
e is difference between two large numbers

 Not important for flow dynamics since p is negligibleBut can cause problems if we want accurate temperatures
since T α eDual Energy Formalism:

grid::ComputePresureDualEnergyFormalism Solution: Also evolve equation for internal energy:

 Select energy to use depending on ratio e/E:

 Select with DualEnergyFormalism = 1
 Use when v/cs > ~20
 Q: Why not just use e?

A: Equation for e is not in conservative form (source term).
 Source term in internal energy equation causes diffusionZeus Solver: grid::ZeusSolver
 Traditional finite difference method

Artificial viscosity (see Stone & Norman 1992)

HydroMethod = 2

 Source step: ZeusSourcePressure (and gravity) update:

Artificial viscosity:Compression heating:Zeus Solver: grid::ZeusSolver
 Transport step: Zeus_xTransporte.g.

Note conservative form (transport part preserves mass) Note vj+1 is face-centered so is really at cell-edge, but density
needs to be interpolated. Zeus uses an upwinded van Leer(linear) interpolation:


Similarly for momentum and energy (and y and z)
 Zeus_yTransport, Zeus_zTransportZeus Solver: grid::ZeusSolver
 PPM is more accurate, slower but Zeus is faster and morerobust.
PPM often fails (“dnu < 0” error) when fast cooling generates
large density gradients.
 Try out new hydro solvers in Enzo 2.0!
 Implementation differences with PPM:Internal energy equation only
 In code, TotalEnergy field is really internal energy (ugh!)Velocities are face-centered

BaryonField[Vel1Num][GINDEX(i,j,k)] really “lives” at i-1/2

II. Block Structured AMRAMR: EvolveHierarchy Root grid NxNxN, so Δx = DomainWidth/N
 Level L defined so Δx = DomainWidth/(N2L)
 Starting with level 0, grid advanced by ΔtMain loop of EvolveHierarchy looks (roughly) like this:

EvolveLevel does the heavy lifting

Time Step: grid::ComputeTimeStep Timestep on level L is minimum of constraints over alllevel L grids:
CourantSafetyFactor

ParticleCourantSafetyFactor

MaximumExpansionFactor

 + others (e.g. MHD, FLD, etc.)

AMR: EvolveLevel
 Levels advanced as follows:Timesteps may not be integer ratios
 (Diagram assumes Courant condition dominates and sound speed

constant so: dt α Δx)This algorithm is defined in EvolveLevelisAdvance grids on level: EvolveLevel
 The logic of EvolveLevel is given (roughly) as:

Already talked about this.

recursiveNext, we’ll talk about these

BC’s: SetBoundaryConditions
 Setting “ghost” zones around outside of domain


grid::SetExternalBoundaryValues

Choices: reflecting, outflow, inflow, periodic
 Only applied to level 0 grids (except periodic)
 Otherwise, two step procedure:Interpolate ghost (boundary) zones from level L-1 grid

grid::InterpolateBoundaryFromParent
 Linear interpolation in time (OldBaryonFields)
 Spatial interpolation controlled by InterpolationMethodSecondOrderA recommended, default (3D, linear in space, monotonic)Copy ghost zones from sibling grids


grid::CheckForOverlap and grid::CopyZonesFromGrid

Projection: grid::ProjectSolutionToParentGrid
 Structured AMR produces redundancy:
coarse and fine grids cover same region

 Need to restore consistency
 Correct coarse cells once grids have all reach the sametime:Flux Correction:

grid::CorrectForRefinedFluxes
 Mismatch of fluxes occurs around boundary of fine grids
Coarse cell just outside boundary used coarse fluxes butcoarse cell inside used fine fluxes Both fine and coarse fluxes saved

from hydro solver

Uncorrectedcoarse value

Coarse flux
across boundary

Sum of fine fluxes
Over 4 (in 3D)
abutting fine cells

Rebuilding the Hierarchy:
RebuildHierarchy Need to check for cells needing more refinement
Check for new grids onlevel 1 (and below)
Check for new grids on
Level 2 (and below)

Check for new grids on
Level 3 (and below)

Refinement Criteria – grid::SetFlaggingField
 Many ways to flag cells for refinement
CellFlaggingMethod = Then rectangular grids must be chosen to cover all

flagged cells with minimum “waste”
Done with machine vision technique Looks for edges (inflection points in number of flagged cells)ProtoSubgrid classIII. Gravity

Self-Gravity (SelfGravity = 1)
 Solve Poisson equationPrepareDensityField


BaryonField[Density] copied to GravitatingMassField



Particle mass is deposited in 8 nearest cells (CIC) Particle position advanced by ½ stepDepositParticleMassField Root grid (level 0):Potential solved with FFTComputePotentialFieldLevelZero

Potential differenced to get acceleration
grid::ComputeAccelerationFieldSelf-Gravity
 Subgrids:Potential interpolated to boundary from parentEach subgrid then solves Poisson equation using multigrid

Grid::PreparePotentialField
Grid::SolveForPotential

Note: this has two issues: Interpolation errors on boundary can propagate to fine levels
 Generally only an issue for steep potentials (point mass)
 Ameliorated by having 6 ghost zones for gravity grid Subgrids can have inconsistent potential gradients across boundary Improved by copying new boundary conditions from sibilings and resolving

the Poisson equation (PotentialIterations = 4 by default) More accurate methods in development

Other Gravitational sources –grid::ComputeAccelerationFieldExternal Can also add fixed potential:UniformGravity – constant fieldPointSourceGravity – single point source

ExternalGravity – NFW profile

IV. Particles

N-body dynamics
 Particles contribute mass to GravitatingMassField Particles accelerated by AccelerationField
Interpolated from grid (from 8 nearest cells) Particles advanced using leapfrog

grid::ComputeAccelerations

 Particles stored in the locally most-refined gridParticlePosition, ParticleVelocity, ParticleMass Tracer particles (massless) also available