-
Notifications
You must be signed in to change notification settings - Fork 0
/
numerical-timestep.tex
150 lines (130 loc) · 6.69 KB
/
numerical-timestep.tex
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
\section{Timestepping}
\label{sec.timestepping}
In \enzo, the integration of the equations being solved is generally
adaptive in time as well as in space. The timestep $\Delta t$ is set
on a level-by-level basis by finding the largest timestep such that
all of the criteria listed below (that are relevant for the simulation
in question) are satisfied. The timestep criteria are given by the
following expressions, showing the one-dimensional case for clarity:
\begin{equation}
\Delta t_{\rm hydro} = \min \left( \kappa_{\rm hydro} \frac{a \Delta x}{c_{s} + |v_x|} \right)_L ,
\label{eqn:dthydro}
\end{equation}
\begin{equation}
\Delta t_{\rm MHD} = \min \left( \kappa_{\rm MHD} \frac{a \Delta x}{v_{f} + |v_x|} \right)_L ,
\label{eqn:dtMHD}
\end{equation}
\begin{equation}
\Delta t_{\rm dm} = \min \left(\kappa_{\rm dm} \frac{a \Delta x}{v_{dm,x}} \right)_L ,
\label{eqn:dtdarkmatter}
\end{equation}
\begin{equation}
\Delta t_{\rm accel} = \min \left( \sqrt{\frac{\Delta x}{|\myvec{g}|}} \right)_L ,
\label{eqn:dtaccel}
\end{equation}
\begin{equation}
\Delta t_{\rm rad} = \min \left( \sqrt{\frac{\Delta x}{|\myvec{a}_{rad}|}} \right)_L,
\label{eqn:dtrad}
\end{equation}
\begin{equation}
\Delta t_{\rm cond} = \min \left( \frac{ k_{\rm cond}}{f_{\rm sp}} \frac{\Delta x^2
n_b}{\kappa_{\rm sp}(T)} \right)_L,
\label{eqn:dtcond}
\end{equation}
\begin{equation}
\Delta t_{\rm exp} = f_{\rm exp} \left( \frac{a}{\dot{a}} \right) ,
\label{eqn:dtexpand}
\end{equation}
In equations~(\ref{eqn:dthydro})-(\ref{eqn:dtcond}), the $\min (
\ldots)_L$ formalism means that this value is calculated for all cells
or particles on a given level L and the minimum overall value is taken
as the timestep.
Equation~(\ref{eqn:dthydro}) ensures that all cells satisfy the
Courant-Freidrichs-Levy (CFL) condition for accuracy and stability of
an explicit finite difference discretization of the Euler equations.
In this equation, $\kappa_{\rm hydro}$ is a dimensionless numerical
constant with a value of $0 < \kappa_{\rm hydro} \leq 1$ (with a
typical value of $\kappa_{\rm hydro} \sim 0.3-0.5$) that ensures that
the CFL condition is always met, and $c_s$ and $v_x$ are the sound
speed and peculiar baryon velocity in a given cell.
Equation~(\ref{eqn:dthydro}) is valid for one dimension, and is used
when the equations of hydrodynamics are being solved. For 2 or 3
dimensions, it was shown by \cite{Godunov1959} that using the harmonic
average of the timestep found along each of the coordinate axes yields
a maximum $\kappa_{\rm hydro} = 0.8$. So letting $\Delta t_x$,
$\Delta t_y$, and $\Delta t_z$ be the analogues of
equation~(\ref{eqn:dthydro}) along the $x,y$ and $z$ axes,
\begin{equation}
\Delta t_{\rm hydro} = \min \left( \frac{\kappa_{\rm hydro}} {1/\Delta t_x
+1/\Delta t_y + 1/\Delta t_z} \right)_L
\end{equation}
For all other criteria except for equation~(\ref{eqn:dtMHD}), multiple
dimensions are accounted for by repeating the one dimensional
criterion along each axis, and taking the minimum.
Equation~(\ref{eqn:dtMHD}) is only enforced when the equations of
magnetohydrodynamics are being solved, and is directly analogous to
equation~(\ref{eqn:dthydro}) in that it ensures that the CFL condition
is being enforced at all times. In this equation, $\kappa_{\rm MHD}$
is a dimensionless numerical constant with a value of $0 < \kappa_{\rm
MHD} \leq 1$ (with a typical value of $\kappa_{\rm MHD} \sim 0.5$)
that ensures that the CFL condition is always met, and $v_f$ and $v_x$
are the ``fast wave speed'' and peculiar baryon velocity in a given
cell. The fast wave speed comes from a stability analysis of the MHD
equations, and is given by:
%
\begin{equation}
v_f = \sqrt{ \frac{1}{2} \left( v_A^2 + c_s^2 + \sqrt{(v_A^2 +
c_s^2)^2 - 4 v_A^2 c_s^2} \right) },
\label{eqn:vfastmhd}
\end{equation}
%
where $c_s$ is the sound speed and $v_A$ is the Alfven speed,
calculated as $v_A = \sqrt{B^2/\rho_B}$ in units where $\mu_0 = 1$.
Equation~(\ref{eqn:dtdarkmatter}) is analogous to
equation~(\ref{eqn:dthydro}) and helps to ensure accuracy in the
N-body solver by requiring that no particle travels more than one cell
width. The parameter $\kappa_{\rm dm}$ is like $\kappa_{\rm hydro}$,
with a similar range of values. This criterion is used when massive
particles are included in a simulation.
Equations~(\ref{eqn:dtaccel}) and~(\ref{eqn:dtrad}) are supplementary
to equation~(\ref{eqn:dthydro}) in that they take into account the
possibility of large accelerations due to either gravity
(equation~\ref{eqn:dtaccel}) or radiation pressure
(equation~\ref{eqn:dtrad}). In equation~(\ref{eqn:dtaccel}),
$\myvec{g}$ is the gravitational acceleration in each cell on level l.
In equation~(\ref{eqn:dtrad}), $\myvec{a}_{\rm rad}$ is the estimated
acceleration due to radiation pressure in each cell on level l,
defined as
\begin{equation}
\myvec{a}_{rad} = \frac{ \sum_i \frac{\dot{E_i}}{c} \hat{r_i} }{m_b}
\end{equation}
where the sum calculates the energy deposited in a cell during the
previous timestep due to \textit{all photon packets} that crossed that
cell, with $\hat{r_i}$ being a unit vector that accounts for the
packet direction.
Equation~(\ref{eqn:dtcond}) is the stability condition for an explicit
solution to the equation of heat conduction. In this expression,
$n_b$ is the baryon number density, $\kappa_{\rm sp}(T)$ is the
Spitzer thermal conductivity, and f$_{\rm sp}$ is a user-defined,
dimensionless conduction suppression factor whose value must be
f$_{\rm sp} \leq 1$. $k_{\rm cond}$ is a dimensionless prefactor
whose value must be $0 < k_{\rm cond} \leq 0.5$, and is exactly 0.5
for the implemention in \enzo. From a practical perspective, it is
useful to note that, unlike other timestep criteria discussed above
(which effectively scale as $\Delta x / \sqrt{T}$ with the grid cell
size $\Delta x$ and temperature $T$), the timestep criterion due to
thermal conduction scales as $\Delta x^2 / T^{2.5}$, which can result
in a rapid decrease in timestep in regions of high resolution and/or
temperature.
Finally, equation~(\ref{eqn:dtexpand}) is a cosmological constraint
that limits the timestep so that the simulated universe only expands
by some fractional amount, f$_{\rm exp}$, during a single step. In
this equation, $a$ and $\dot{a}$ refer to the scale factor of the
universe and its rate of change, respectively. This criteria is
necessary because the expansion of the universe and its first
derivative with respect to time both appear in the equations of
cosmological (magneto)hydrodynamics and particle motion, and some
limit is required for the stability of the PPM algorithm in comoving
coordinates. This criterion typically limits the simulation's
timestep only during the earliest phases of a cosmological simulation,
before substantial structure has formed.