Commit 28265917 authored by Matthieu Schaller's avatar Matthieu Schaller

Merge branch 'cosmological_integration' into 'master'

Cosmological integration

Closes #405

See merge request !979
parents 856bcb95 583944b3
......@@ -113,6 +113,7 @@ tests/testSelectOutput
tests/testReading
tests/testSingle
tests/testTimeIntegration
tests/testIntegration
tests/testSPHStep
tests/testKernel
tests/testKernelGrav
......
This diff is collapsed.
......@@ -163,6 +163,15 @@ struct cosmology {
/*! Log of final expansion factor */
double log_a_end;
/*! Log of starting expansion factor of the interpolation tables */
double log_a_table_begin;
/*! Log of final expansion factor of the interpolation tables */
double log_a_table_end;
/*! Scale-factor interpolation table */
double *log_a_interp_table;
/*! Drift factor interpolation table */
double *drift_fac_interp_table;
......@@ -175,15 +184,9 @@ struct cosmology {
/*! Kick factor (hydro correction) interpolation table (GIZMO-MFV only) */
double *hydro_kick_corr_interp_table;
/*! Time interpolation table */
/*! Time since Big Bang interpolation table */
double *time_interp_table;
/*! Scale factor interpolation table */
double *scale_factor_interp_table;
/*! Time between Big Bang and first entry in the table */
double time_interp_table_offset;
/*! Time at the present-day (a=1) */
double universe_age_at_present_day;
};
......@@ -191,6 +194,8 @@ struct cosmology {
void cosmology_update(struct cosmology *c, const struct phys_const *phys_const,
integertime_t ti_current);
double cosmology_get_scale_factor(const struct cosmology *cosmo,
integertime_t ti);
double cosmology_get_drift_factor(const struct cosmology *cosmo,
const integertime_t ti_start,
const integertime_t ti_end);
......@@ -214,7 +219,8 @@ double cosmology_get_delta_time_from_scale_factors(const struct cosmology *c,
const double a_start,
const double a_end);
double cosmology_get_scale_factor(const struct cosmology *cosmo, double t);
double cosmology_get_scale_factor_from_time(const struct cosmology *cosmo,
double t);
double cosmology_get_time_since_big_bang(const struct cosmology *c, double a);
void cosmology_init(struct swift_params *params, const struct unit_system *us,
......
......@@ -106,7 +106,7 @@ void output_list_read_file(struct output_list *outputlist, const char *filename,
if (type == OUTPUT_LIST_REDSHIFT) *time = 1. / (1. + *time);
if (cosmo && type == OUTPUT_LIST_AGE)
*time = cosmology_get_scale_factor(cosmo, *time);
*time = cosmology_get_scale_factor_from_time(cosmo, *time);
/* Update size */
ind += 1;
......
......@@ -28,7 +28,7 @@ TESTS = testGreetings testMaths testReading.sh testKernel \
testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf \
testPotentialPair testEOS testUtilities testSelectOutput.sh \
testCbrt testCosmology testOutputList testFormat.sh \
testCbrt testCosmology testIntegration testOutputList testFormat.sh \
test27cellsStars.sh test27cellsStarsPerturbed.sh testHydroMPIrules
# List of test programs to compile
......@@ -40,7 +40,7 @@ check_PROGRAMS = testGreetings testReading testTimeIntegration \
testRiemannHLLC testMatrixInversion testDump testLogger \
testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
testGravityDerivatives testPotentialSelf testPotentialPair testEOS testUtilities \
testSelectOutput testCbrt testCosmology testOutputList test27cellsStars \
testSelectOutput testCbrt testIntegration testCosmology testOutputList test27cellsStars \
test27cellsStars_subset testCooling testComovingCooling testFeedback testHashmap testHydroMPIrules
# Rebuild tests when SWIFT is updated.
......@@ -59,6 +59,8 @@ testReading_SOURCES = testReading.c
testSelectOutput_SOURCES = testSelectOutput.c
testIntegration_SOURCES = testIntegration.c
testCosmology_SOURCES = testCosmology.c
testOutputList_SOURCES = testOutputList.c
......
......@@ -20,11 +20,15 @@
/* Some standard headers. */
#include "../config.h"
/* Some standard headers */
#include <fenv.h>
#include <math.h>
/* Includes. */
#include "swift.h"
#define N_CHECK 20
#define TOLERANCE 1e-7
#define N_CHECK 200
#define TOLERANCE 1e-6
void test_params_init(struct swift_params *params) {
parser_init("", params);
......@@ -32,12 +36,21 @@ void test_params_init(struct swift_params *params) {
parser_set_param(params, "Cosmology:Omega_lambda:0.6910");
parser_set_param(params, "Cosmology:Omega_b:0.0486");
parser_set_param(params, "Cosmology:h:0.6774");
parser_set_param(params, "Cosmology:a_begin:0.1");
parser_set_param(params, "Cosmology:a_begin:0.01");
parser_set_param(params, "Cosmology:a_end:1.0");
}
int main(int argc, char *argv[]) {
/* Initialize CPU frequency, this also starts time. */
unsigned long long cpufreq = 0;
clocks_set_cpufreq(cpufreq);
/* Choke on FP-exceptions */
#ifdef HAVE_FE_ENABLE_EXCEPT
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
#endif
message("Initialization...");
/* pseudo initialization of params */
......@@ -59,15 +72,16 @@ int main(int argc, char *argv[]) {
message("Start checking computation...");
for (int i = 0; i < N_CHECK; i++) {
double a = 0.1 + 0.9 * i / (N_CHECK - 1.);
double a = 0.01 + 0.99 * i / (N_CHECK - 1.);
/* Compute a(t(a)) and check if same results */
double tmp = cosmology_get_time_since_big_bang(&cosmo, a);
tmp = cosmology_get_scale_factor(&cosmo, tmp);
double time = cosmology_get_time_since_big_bang(&cosmo, a);
double my_a = cosmology_get_scale_factor_from_time(&cosmo, time);
/* check accuracy */
tmp = (tmp - a) / a;
message("Accuracy of %g at a=%g", tmp, a);
assert(fabs(tmp) < TOLERANCE);
double rel_err = (my_a - a) / a;
message("Accuracy of %g at a=%g", rel_err, a);
assert(fabs(rel_err) < TOLERANCE);
}
message("Everything seems fine with cosmology.");
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2019 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#include "../config.h"
/* Some standard headers */
#include <fenv.h>
#include <math.h>
/* Local headers. */
#include "swift.h"
#ifdef HAVE_LIBGSL
#include <gsl/gsl_integration.h>
#endif
static INLINE double w_tilde(double a, double w0, double wa) {
return (a - 1.) * wa - (1. + w0 + wa) * log(a);
}
static INLINE double E(double Or, double Om, double Ok, double Ol, double w0,
double wa, double a) {
const double a_inv = 1. / a;
return sqrt(Or * a_inv * a_inv * a_inv * a_inv + Om * a_inv * a_inv * a_inv +
Ok * a_inv * a_inv + Ol * exp(3. * w_tilde(a, w0, wa)));
}
static INLINE double drift_integrand(double a, void *param) {
const struct cosmology *c = (const struct cosmology *)param;
const double Omega_r = c->Omega_r;
const double Omega_m = c->Omega_m;
const double Omega_k = c->Omega_k;
const double Omega_l = c->Omega_lambda;
const double w_0 = c->w_0;
const double w_a = c->w_a;
const double H0 = c->H0;
const double a_inv = 1. / a;
const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w_0, w_a, a);
const double H = H0 * E_z;
return (1. / H) * a_inv * a_inv * a_inv;
}
int main(int argc, char *argv[]) {
/* Initialize CPU frequency, this also starts time. */
unsigned long long cpufreq = 0;
clocks_set_cpufreq(cpufreq);
/* Choke on FP-exceptions */
#ifdef HAVE_FE_ENABLE_EXCEPT
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
#endif
/* Fake parameter file with Planck+13 cosmology
* and the usual cosmological system of units */
struct swift_params *params = malloc(sizeof(struct swift_params));
parser_set_param(params, "InternalUnitSystem:UnitMass_in_cgs:1.98848e43");
parser_set_param(params,
"InternalUnitSystem:UnitLength_in_cgs:3.08567758e24");
parser_set_param(params, "InternalUnitSystem:UnitVelocity_in_cgs:1e5");
parser_set_param(params, "InternalUnitSystem:UnitCurrent_in_cgs:1.");
parser_set_param(params, "InternalUnitSystem:UnitTemp_in_cgs:1.");
parser_set_param(params, "Cosmology:Omega_m:0.307");
parser_set_param(params, "Cosmology:Omega_lambda:0.693");
parser_set_param(params, "Cosmology:Omega_b:0.0482519");
parser_set_param(params, "Cosmology:h:0.6777");
parser_set_param(params, "Cosmology:a_begin:0.0078125");
parser_set_param(params, "Cosmology:a_end:1.");
/* Initialise everything */
struct unit_system units;
units_init_from_params(&units, params, "InternalUnitSystem");
struct phys_const phys_const;
phys_const_init(&units, params, &phys_const);
struct cosmology cosmo;
cosmology_init(params, &units, &phys_const, &cosmo);
/* Initalise the GSL workspace */
size_t workspace_size = 100000;
gsl_integration_workspace *workspace =
gsl_integration_workspace_alloc(workspace_size);
gsl_function F = {&drift_integrand, &cosmo};
/* Loop over all the reasonable time-bins */
for (int bin = 6; bin < 24; ++bin) {
const int num_steps = (1LL << bin);
const integertime_t time_step_size = max_nr_timesteps / num_steps;
double min_err = 0.;
double max_err = 0.;
double sum_err = 0.;
integertime_t ti_min = 0;
integertime_t ti_max = 0;
message("Testing %d steps", num_steps);
/* Cycle through the time-steps */
for (integertime_t ti = 0; ti < max_nr_timesteps; ti += time_step_size) {
const integertime_t ti_beg = ti;
const integertime_t ti_end =
min(ti + time_step_size, max_nr_timesteps - 1);
const double a_beg = cosmology_get_scale_factor(&cosmo, ti_beg);
const double a_end = cosmology_get_scale_factor(&cosmo, ti_end);
/* Get the drift factor from SWIFT */
const double swift_drift_fac =
cosmology_get_drift_factor(&cosmo, ti, ti + time_step_size);
/* Get the exact drift factor */
double exact_drift_fac = 0., abserr;
gsl_integration_qag(&F, a_beg, a_end, 0, 1.0e-12, workspace_size,
GSL_INTEG_GAUSS61, workspace, &exact_drift_fac,
&abserr);
const double rel_err = 0.5 * (swift_drift_fac - exact_drift_fac) /
(swift_drift_fac + exact_drift_fac);
if (rel_err > max_err) {
max_err = rel_err;
ti_max = ti;
}
if (rel_err < min_err) {
min_err = rel_err;
ti_min = ti;
}
sum_err += fabs(rel_err);
}
message("Max error: %14e at a=[%.9f %.9f] ", max_err,
cosmology_get_scale_factor(&cosmo, ti_max),
cosmology_get_scale_factor(&cosmo, ti_max + time_step_size));
message("Min error: %14e at a=[%.9f %.9f]", min_err,
cosmology_get_scale_factor(&cosmo, ti_min),
cosmology_get_scale_factor(&cosmo, ti_min + time_step_size));
message("Sum error: %14e", sum_err);
message("Mean error: %14e", sum_err / num_steps);
if (max_err > 1e-8 || min_err < -1e-8)
error("Error too large to be acceptable");
}
#ifdef HAVE_LIBGSL
return 0;
#else
return 0;
#endif
}
......@@ -19,6 +19,10 @@
******************************************************************************/
#include "../config.h"
/* Some standard headers */
#include <fenv.h>
#include <math.h>
/* Includes. */
#include "swift.h"
......@@ -86,7 +90,8 @@ void test_cosmo(struct engine *e, char *name, int with_assert) {
double output_time = 0;
/* Test Time */
e->time_base = log(e->time_end / e->cosmology->a_begin) / max_nr_timesteps;
e->time_base =
log(e->cosmology->a_end / e->cosmology->a_begin) / max_nr_timesteps;
e->ti_current = 0;
e->policy = engine_policy_cosmology;
......@@ -115,6 +120,16 @@ void test_cosmo(struct engine *e, char *name, int with_assert) {
};
int main(int argc, char *argv[]) {
/* Initialize CPU frequency, this also starts time. */
unsigned long long cpufreq = 0;
clocks_set_cpufreq(cpufreq);
/* Choke on FP-exceptions */
#ifdef HAVE_FE_ENABLE_EXCEPT
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
#endif
/* Create a structure to read file into. */
struct swift_params params;
......
\subsection{Choice of co-moving coordinates}
\label{ssec:ccordinates}
Note that, unlike the \gadget convention, we do not express quantities with
``little h'' ($h$) included; for instance units of length are expressed in
units of $\rm{Mpc}$ and not ${\rm{Mpc}}/h$. Similarly, the time integration
operators (see below) also include an $h$-factor via the explicit appearance of
Note that, unlike the convention used by many codes, we do not express
quantities with ``little h'' ($h$) included; for instance units of
length are expressed in units of $\rm{Mpc}$ and not
${\rm{Mpc}}/h$. Similarly, the time integration operators (see
below)include an $h$-factor via the explicit appearance of
the Hubble constant.\\
In physical coordinates, the Lagrangian for a particle $i$ in the
\cite{Springel2002} flavour of SPH with gravity reads
\begin{equation}
......@@ -61,82 +63,8 @@ codes\footnote{One additional inconvenience of our choice of
scale-factor. The signal velocity entering the time-step calculation
will hence read
$v_{\rm sig} = a\dot{\mathbf{r}'} + c = \frac{1}{a} \left(
|\mathbf{v}'| + a^{(5 - 3\gamma)/2}c'\right)$.}.
This choice implies that $\dot{\mathbf{v}'} = 2H\mathbf{v}' + a^2\ddot{\mathbf{r}'}$.
\subsubsection{SPH equations}
Using the SPH definition of density,
$\rho_i' = \sum_jm_jW(\mathbf{r}_{j}'-\mathbf{r}_{i}',h_i') =
\sum_jm_jW_{ij}'(h_i')$, we can follow \cite{Price2012} and apply the
Euler-Lagrange equations to write
\begin{alignat}{3}
\dot{\mathbf{r}}_i'&= \frac{1}{a^2} \mathbf{v}_i'& \label{eq:cosmo_eom_r} \\
\dot{\mathbf{v}}_i' &= -\sum_j m_j &&\left[\frac{1}{a^{3(\gamma-1)}}f_i'A_i'\rho_i'^{\gamma-2}\mathbf{\nabla}_i'W_{ij}'(h_i)\right. \nonumber\\
& && + \left. \frac{1}{a^{3(\gamma-1)}}f_j'A_j'\rho_j'^{\gamma-2}\mathbf{\nabla}_i'W_{ij}'(h_j)\right. \nonumber\\
& && + \left. \frac{1}{a}\mathbf{\nabla}_i'\phi'\right] \label{eq:cosmo_eom_v}
\end{alignat}
with
\begin{equation}
f_i' = \left[1 + \frac{h_i'}{3\rho_i'}\frac{\partial
\rho_i'}{\partial h_i'}\right]^{-1}, \qquad \mathbf{\nabla}_i'
\equiv \frac{\partial}{\partial \mathbf{r}_{i}'}. \nonumber
\end{equation}
These correspond to the equations of motion for density-entropy SPH
\citep[e.g. eq. 14 of][]{Hopkins2013} with cosmological and
gravitational terms. SPH flavours that evolve the internal energy $u$ instead of the
entropy require the additional equation of motion describing the evolution of
$u'$:
\begin{equation}
\dot{u}_i' = \frac{1}{a^2}\frac{P_i'}{\rho_i'^2} f_i'\sum_jm_j\left(\mathbf{v}_i' -
\mathbf{v}_j'\right)\cdot\mathbf{\nabla}_i'W_{ij}'(h_i).
\label{eq:cosmo_eom_u}
\end{equation}
|\mathbf{v}'| + a^{(5 - 3\gamma)/2}c'\right)$.}. For the
derivatives, this choice implies
$\dot{\mathbf{v}'} = 2H\mathbf{v}' + a^2\ddot{\mathbf{r}'}$.
In all these cases, the scale-factors appearing in the equations are
later absorbed in the time-integration operators
(Sec.~\ref{ssec:operators}) such that the RHS of the equations of
motions is identical for the primed quantities to the ones obtained in
the non-cosmological case for the physical quantities.
Additional terms in the SPH equations of motion (e.g. viscosity
switches) often rely on the velocity divergence and curl. Their SPH
estimators $\langle\cdot\rangle$ in physical coordinates can be
related to their estimators based on our primed-coordinates using:
\begin{align}
\left\langle \mathbf{\nabla}\cdot\dot{\mathbf{r}}_i \right\rangle &=
\frac{1}{a^2} \left\langle
\mathbf{\nabla}'\cdot\mathbf{v}_i'\right\rangle =
\frac{1}{a^2\rho_i'}\sum_j m_j\left(\mathbf{v}_j' -
\mathbf{v}_i'\right) \cdot \mathbf{\nabla}_i'W_{ij}'(h_i), \nonumber \\
\left\langle \mathbf{\nabla}\times\dot{\mathbf{r}}_i \right\rangle &=
\frac{1}{a^2} \left\langle
\mathbf{\nabla}'\times\mathbf{v}_i'\right\rangle =
\frac{1}{a^2\rho_i'}\sum_j m_j\left(\mathbf{v}_j' -
\mathbf{v}_i'\right) \times \mathbf{\nabla}_i'W_{ij}'(h_i). \nonumber
\end{align}
We finally give the expressions for the \cite{Monaghan1997} viscosity
term, that enters most SPH flavours, in our system of coordinates. The
viscosity ``tensor'' $\Pi_{ij} \equiv \frac{1}{2}\alpha_{\rm visc} v_{\rm
sig}'\mu_{ij}/\left(\rho_i + \rho_j\right)$ can be computed
in the primmed coordinates from the following quantities
\begin{align}
\omega_{ij}' &= \left(\mathbf{v}_i' - \mathbf{v}_j'\right) \cdot
\left(\mathbf{r}_i' - \mathbf{r}_j'\right), \\
\mu_{ij}' &=
a^{(3\gamma-5)/2} \left(\omega_{ij}' + a^2H\left|\mathbf{r}_i' -
\mathbf{r}_j'\right|^2 \right) / |\mathbf{r}_i' - \mathbf{r}_j' |,
\\
v_{\rm sig}' &= c_i' + c_j' - 3 \mu_{ij}'.
\end{align}
which leads to $\Pi_{ij}'=a^{3\gamma}\Pi_{ij}$. Note that he last quantity is
also used to compute the time-step size. The evolution of entropy is
then
\begin{equation}
\dot{A}_i = \dot{A}_i' = \frac{1}{a^2}\frac{1}{2}\frac{\gamma-1}{\rho_i'^{\gamma-1}} \sum_j
m_j \Pi_{ij}' \left(\mathbf{v}_i' -
\mathbf{v}_j'\right)\cdot\mathbf{\nabla}_i'W_{ij}'(h_i),
\end{equation}
indicating that the entropy evolves with the same scale-factor
dependence as the comoving positions (eq.~\ref{eq:cosmo_eom_r}).
......@@ -26,7 +26,7 @@
\label{firstpage}
\begin{abstract}
Making cosmology great again.
Notes on the cosmology equations evolved in \swift
\end{abstract}
\begin{keywords}
......@@ -39,13 +39,15 @@ Making cosmology great again.
\input{coordinates}
\input{sph}
\input{operators}
\input{timesteps}
\input{artificialvisc}
%\input{artificialvisc}
\input{gizmo}
%\input{gizmo}
\bibliographystyle{mnras}
\bibliography{./bibliography.bib}
......
......@@ -48,7 +48,7 @@ of $a_{\rm age}$ equally spaced between $\log(a_{\rm start})$ and $\log(a_{\rm
end})$. The values are obtained via adaptive quadrature using the 61-points
Gauss-Konrod rule implemented in the {\sc gsl} library \citep{GSL} with a
relative error limit of $\epsilon=10^{-10}$. The value for a specific $a$ (over
the course of a simulation run) is then obtained by linear interpolation of the
the course of a simulation run) is then obtained by cubic interpolation of the
table.
\subsubsection{Additional quantities}
......
......@@ -47,11 +47,12 @@ and the corresponding change in redshift:
\Delta z &= \frac{1}{a_n} - \frac{1}{a_{n+1}} \approx -\frac{H}{a} \Delta t_{\rm cosmic}.
\end{align}
Following the same method as for the age of the Universe
(sec. \ref{ssec:flrw}), these three non-trivial integrals are
evaluated numerically at the start of the simulation for a series
$10^4$ values of $a$ placed at regular intervals between $\log a_{\rm
begin}$ and $\log a_{\rm end}$. The values for a specific pair of
scale-factors $a_n$ and $a_{n+1}$ are then obtained by interpolating
that table linearly.
(sec. \ref{ssec:flrw}), these three non-trivial integrals are evaluated
numerically at the start of the simulation for a series $10^4$ values of
$a$ placed at regular intervals between $\log a_{\rm begin}$ and
$\log a_{\rm end}$. The values for a specific pair of scale-factors $a_n$
and $a_{n+1}$ are then obtained by interpolating that table with a cubic
polynomial. An accuracy better than one part in $10^8$ is achieved for all
reasonable cosmologies and time-step lengths.
\subsubsection{SPH equations}
Using the SPH definition of density,
$\rho_i' = \sum_jm_jW(\mathbf{r}_{j}'-\mathbf{r}_{i}',h_i') =
\sum_jm_jW_{ij}'(h_i')$, we can follow \cite{Price2012} and apply the
Euler-Lagrange equations to write
\begin{alignat}{3}
\dot{\mathbf{r}}_i'&= \frac{1}{a^2} \mathbf{v}_i'& \label{eq:cosmo_eom_r} \\
\dot{\mathbf{v}}_i' &= -\sum_j m_j &&\left[\frac{1}{a^{3(\gamma-1)}}f_i'A_i'\rho_i'^{\gamma-2}\mathbf{\nabla}_i'W_{ij}'(h_i)\right. \nonumber\\
& && + \left. \frac{1}{a^{3(\gamma-1)}}f_j'A_j'\rho_j'^{\gamma-2}\mathbf{\nabla}_i'W_{ij}'(h_j)\right. \nonumber\\
& && + \left. \frac{1}{a}\mathbf{\nabla}_i'\phi'\right] \label{eq:cosmo_eom_v}
\end{alignat}
with
\begin{equation}
f_i' = \left[1 + \frac{h_i'}{3\rho_i'}\frac{\partial
\rho_i'}{\partial h_i'}\right]^{-1}, \qquad \mathbf{\nabla}_i'
\equiv \frac{\partial}{\partial \mathbf{r}_{i}'}. \nonumber
\end{equation}
These correspond to the equations of motion for density-entropy SPH
\citep[e.g. eq. 14 of][]{Hopkins2013} with cosmological and
gravitational terms. SPH flavours that evolve the internal energy $u$ instead of the
entropy require the additional equation of motion describing the evolution of
$u'$:
\begin{equation}
\dot{u}_i' = \frac{1}{a^2}\frac{P_i'}{\rho_i'^2} f_i'\sum_jm_j\left(\mathbf{v}_i' -
\mathbf{v}_j'\right)\cdot\mathbf{\nabla}_i'W_{ij}'(h_i).
\label{eq:cosmo_eom_u}
\end{equation}
In all these cases, the scale-factors appearing as the first term on
the RHS of the equations are later absorbed in the time-integration
operators (Sec.~\ref{ssec:operators}) such that the RHS of the
equations of motions is identical for the primed quantities to the
ones obtained in the non-cosmological case for the physical
quantities.
Additional terms in the SPH equations of motion (e.g. viscosity
switches) often rely on the velocity divergence and curl. Their SPH
estimators $\langle\cdot\rangle$ in physical coordinates can be
related to their estimators based on our primed-coordinates using:
\begin{align}
\left\langle \mathbf{\nabla}\cdot\dot{\mathbf{r}}_i \right\rangle &=
\frac{1}{a^2} \left\langle
\mathbf{\nabla}'\cdot\mathbf{v}_i'\right\rangle =
\frac{1}{a^2\rho_i'}\sum_j m_j\left(\mathbf{v}_j' -
\mathbf{v}_i'\right) \cdot \mathbf{\nabla}_i'W_{ij}'(h_i), \nonumber \\
\left\langle \mathbf{\nabla}\times\dot{\mathbf{r}}_i \right\rangle &=
\frac{1}{a^2} \left\langle
\mathbf{\nabla}'\times\mathbf{v}_i'\right\rangle =
\frac{1}{a^2\rho_i'}\sum_j m_j\left(\mathbf{v}_j' -
\mathbf{v}_i'\right) \times \mathbf{\nabla}_i'W_{ij}'(h_i). \nonumber
\end{align}
We finally give the expressions for the \cite{Monaghan1997} viscosity
term, that enters most SPH flavours, in our system of coordinates. The
viscosity ``tensor'' $\Pi_{ij} \equiv \frac{1}{2}\alpha_{\rm visc} v_{\rm
sig}'\mu_{ij}/\left(\rho_i + \rho_j\right)$ can be computed
in the primmed coordinates from the following quantities
\begin{align}
\omega_{ij}' &= \left(\mathbf{v}_i' - \mathbf{v}_j'\right) \cdot
\left(\mathbf{r}_i' - \mathbf{r}_j'\right), \\
\mu_{ij}' &=
a^{(3\gamma-5)/2} \left(\omega_{ij}' + a^2H\left|\mathbf{r}_i' -
\mathbf{r}_j'\right|^2 \right) / |\mathbf{r}_i' - \mathbf{r}_j' |,
\\
v_{\rm sig}' &= c_i' + c_j' - 3 \mu_{ij}'.
\end{align}
which leads to $\Pi_{ij}'=a^{3\gamma}\Pi_{ij}$. Note that he last quantity is
also used to compute the time-step size. The evolution of entropy is
then
\begin{equation}
\dot{A}_i = \dot{A}_i' = \frac{1}{a^2}\frac{1}{2}\frac{\gamma-1}{\rho_i'^{\gamma-1}} \sum_j
m_j \Pi_{ij}' \left(\mathbf{v}_i' -
\mathbf{v}_j'\right)\cdot\mathbf{\nabla}_i'W_{ij}'(h_i),
\end{equation}
indicating that the entropy evolves with the same scale-factor
dependence as the comoving positions (eq.~\ref{eq:cosmo_eom_r}).
......@@ -44,3 +44,5 @@ dominates the overall time-step size calculation.
\int_{a_n}^{a_{n+1}} H dt = \int_{a_n}^{a_{n+1}} \frac{da}{a} =
\log{a_{n+1}} - \log{a_n}.
\end{equation}
\subsubsection{Maximal time-step size}
\ No newline at end of file
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment