Commit 51642837 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added cosmological terms to the Gadget2 hydro scheme.

parent 27623509
......@@ -2392,7 +2392,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
} else if (!c->split && force && ti_current > ti_old_part) {
/* Drift from the last time the cell was drifted to the current time */
double dt_drift, dt_kick_grav, dt_kick_hydro;
double dt_drift, dt_kick_grav, dt_kick_hydro, dt_therm;
if (e->policy & engine_policy_cosmology) {
dt_drift =
cosmology_get_drift_factor(e->cosmology, ti_old_part, ti_current);
......@@ -2400,10 +2400,13 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
cosmology_get_grav_kick_factor(e->cosmology, ti_old_part, ti_current);
dt_kick_hydro = cosmology_get_hydro_kick_factor(e->cosmology, ti_old_part,
ti_current);
dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_old_part,
ti_current);
} else {
dt_drift = (ti_current - ti_old_part) * e->time_base;
dt_kick_grav = (ti_current - ti_old_part) * e->time_base;
dt_kick_hydro = (ti_current - ti_old_part) * e->time_base;
dt_therm = (ti_current - ti_old_part) * e->time_base;
}
/* Loop over all the gas particles in the cell */
......@@ -2415,8 +2418,8 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
struct xpart *const xp = &xparts[k];
/* Drift... */
drift_part(p, xp, dt_drift, dt_kick_hydro, dt_kick_grav, ti_old_part,
ti_current);
drift_part(p, xp, dt_drift, dt_kick_hydro, dt_kick_grav, dt_therm,
ti_old_part, ti_current);
/* Limit h to within the allowed range */
p->h = min(p->h, hydro_h_max);
......
......@@ -136,12 +136,15 @@ void cosmology_update(struct cosmology *c, integertime_t ti_current) {
const double a_inv = 1. / a;
c->a = a;
c->a_inv = a_inv;
c->a2_inv = a_inv * a_inv;
c->a3_inv = a_inv * a_inv * a_inv;
c->a_factor_sound_speed =
pow(a, -1.5 * hydro_gamma_minus_one); /* a^{3*(1-gamma)/2} */
c->a_factor_grav_accel = a_inv * a_inv; /* 1 / a^2 */
c->a_factor_hydro_accel =
powf(a, -3.f * hydro_gamma + 2.f); /* 1 / a^(3*gamma - 2) */
pow(a, -3. * hydro_gamma + 2.); /* 1 / a^(3*gamma - 2) */
c->a_factor_mu =
pow(a, 0.5 * (3. * hydro_gamma - 5.)); /* a^{(3*gamma - 5) / 2} */
/* Redshift */
c->z = a_inv - 1.;
......@@ -457,8 +460,10 @@ void cosmology_init_no_cosmo(const struct swift_params *params,
c->H = 1.;
c->a = 1.;
c->a_inv = 1.;
c->a2_inv = 1.;
c->a3_inv = 1.;
c->a_factor_sound_speed = 1.;
c->a_factor_mu = 1.;
c->a_factor_hydro_accel = 1.;
c->a_factor_grav_accel = 1.;
......
......@@ -38,12 +38,18 @@ struct cosmology {
/*! Inverse of the current expansion factor of the Universe */
double a_inv;
/*! Inverse square of the current expansion factor of the Universe */
double a2_inv;
/*! Inverse cube of the current expansion factor of the Universe */
double a3_inv;
/*! Power of the scale-factor used for sound-speed conversion */
double a_factor_sound_speed;
/*! Power of the scale-factor used for relative velocities in viscosity ter */
double a_factor_mu;
/*! Power of the scale-factor used for gravity accelerations */
double a_factor_grav_accel;
......
......@@ -71,13 +71,14 @@ __attribute__((always_inline)) INLINE static void drift_gpart(
* @param dt_drift The drift time-step
* @param dt_kick_grav The kick time-step for gravity accelerations.
* @param dt_kick_hydro The kick time-step for hydro accelerations.
* @param dt_therm The drift time-step for thermodynamic quantities.
* @param ti_old Integer start of time-step (for debugging checks).
* @param ti_current Integer end of time-step (for debugging checks).
*/
__attribute__((always_inline)) INLINE static void drift_part(
struct part *restrict p, struct xpart *restrict xp, double dt_drift,
double dt_kick_hydro, double dt_kick_grav, integertime_t ti_old,
integertime_t ti_current) {
double dt_kick_hydro, double dt_kick_grav, double dt_therm,
integertime_t ti_old, integertime_t ti_current) {
#ifdef SWIFT_DEBUG_CHECKS
if (p->ti_drift != ti_old)
......@@ -104,7 +105,7 @@ __attribute__((always_inline)) INLINE static void drift_part(
p->v[2] += xp->a_grav[2] * dt_kick_grav;
/* Predict the values of the extra fields */
hydro_predict_extra(p, xp, dt_drift);
hydro_predict_extra(p, xp, dt_drift, dt_therm);
/* Compute offsets since last cell construction */
for (int k = 0; k < 3; k++) {
......
......@@ -282,9 +282,10 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* and add the self-contribution term.
*
* @param p The particle to act upon
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
struct part *restrict p) {
struct part *restrict p, const struct cosmology *cosmo) {
/* Some smoothing length multiples. */
const float h = p->h;
......@@ -305,14 +306,15 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.wcount_dh *= h_inv_dim_plus_one;
const float rho_inv = 1.f / p->rho;
const float a_inv2 = cosmo->a_inv * cosmo->a_inv;
/* Finish calculation of the velocity curl components */
p->density.rot_v[0] *= h_inv_dim_plus_one * rho_inv;
p->density.rot_v[1] *= h_inv_dim_plus_one * rho_inv;
p->density.rot_v[2] *= h_inv_dim_plus_one * rho_inv;
/* Finish calculation of the (physical) velocity curl components */
p->density.rot_v[0] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
p->density.rot_v[1] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
/* Finish calculation of the velocity divergence */
p->density.div_v *= h_inv_dim_plus_one * rho_inv;
/* Finish calculation of the (physical) velocity divergence */
p->density.div_v *= h_inv_dim_plus_one * a_inv2 * rho_inv;
}
/**
......@@ -320,9 +322,11 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
struct part *restrict p, struct xpart *restrict xp) {
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) {
/* Some smoothing length multiples. */
const float h = p->h;
......@@ -347,13 +351,13 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
* @param ti_current The current time (on the timeline)
* @param timeBase The minimal time-step size
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp) {
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) {
const float fac_mu = 1.f; /* Will change with cosmological integration */
const float fac_mu = cosmo->a_factor_mu;
/* Inverse of the physical density */
const float rho_inv = 1.f / p->rho;
......@@ -377,7 +381,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the Balsara switch */
const float balsara =
abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);
abs_div_v / (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
/* Compute the "grad h" term */
const float omega_inv =
......@@ -443,18 +447,17 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
*
* @param p The particle
* @param xp The extended data of the particle
* @param dt The drift time-step.
* @param t0 The time at the start of the drift
* @param t1 The time at the end of the drift
* @param timeBase The minimal time-step size
* @param dt_drift The drift time-step for positions.
* @param dt_therm The drift time-step for thermal quantities.
*/
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part *restrict p, const struct xpart *restrict xp, float dt) {
struct part *restrict p, const struct xpart *restrict xp, float dt_drift,
float dt_therm) {
const float h_inv = 1.f / p->h;
/* Predict smoothing length */
const float w1 = p->force.h_dt * h_inv * dt;
const float w1 = p->force.h_dt * h_inv * dt_drift;
if (fabsf(w1) < 0.2f)
p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
else
......@@ -468,7 +471,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->rho *= expf(w2);
/* Predict the entropy */
p->entropy += p->entropy_dt * dt;
p->entropy += p->entropy_dt * dt_therm;
/* Re-compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
......@@ -491,14 +494,15 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
* Multiplies the forces and accelerationsby the appropiate constants
*
* @param p The particle to act upon
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part *restrict p) {
struct part *restrict p, const struct cosmology *cosmo) {
p->force.h_dt *= p->h * hydro_dimension_inv;
p->entropy_dt =
0.5f * gas_entropy_from_internal_energy(p->rho, p->entropy_dt);
p->entropy_dt = 0.5f * cosmo->a2_inv *
gas_entropy_from_internal_energy(p->rho, p->entropy_dt);
}
/**
......
......@@ -406,7 +406,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
float wi, wj, wi_dx, wj_dx;
const float fac_mu = 1.f; /* Will change with cosmological integration */
/* Will change with cosmological integration */
const float fac_mu = 1.f;
const float a2_Hubble = 0.f;
const float r = sqrtf(r2);
const float r_inv = 1.0f / r;
......@@ -446,7 +448,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
(pi->v[2] - pj->v[2]) * dx[2];
(pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
/* Balsara term */
const float balsara_i = pi->force.balsara;
......@@ -513,7 +515,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
float wi, wj, wi_dx, wj_dx;
const float fac_mu = 1.f; /* Will change with cosmological integration */
/* Will change with cosmological integration */
const float fac_mu = 1.f;
const float a2_Hubble = 0.f;
const float r = sqrtf(r2);
const float r_inv = 1.0f / r;
......@@ -553,7 +557,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
(pi->v[2] - pj->v[2]) * dx[2];
(pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
/* Balsara term */
const float balsara_i = pi->force.balsara;
......
......@@ -653,6 +653,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
struct xpart *restrict xparts = c->xparts;
const struct engine *e = r->e;
const struct space *s = e->s;
const struct hydro_space *hs = &s->hs;
const struct cosmology *cosmo = e->cosmology;
const float hydro_h_max = e->hydro_properties->h_max;
const float eps = e->hydro_properties->h_tolerance;
const float hydro_eta_dim =
......@@ -713,7 +715,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
} else {
/* Finish the density calculation */
hydro_end_density(p);
hydro_end_density(p, cosmo);
/* Compute one step of the Newton-Raphson scheme */
const float n_sum = p->density.wcount * h_old_dim;
......@@ -750,10 +752,11 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
redo += 1;
/* Re-initialise everything */
hydro_init_part(p, &s->hs);
hydro_init_part(p, hs);
/* Off we go ! */
continue;
} else {
/* Ok, this particle is a lost cause... */
......@@ -761,7 +764,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* Do some damage control if no neighbours at all were found */
if (p->density.wcount == kernel_root * kernel_norm)
hydro_part_has_no_neighbours(p, xp);
hydro_part_has_no_neighbours(p, xp, cosmo);
}
}
......@@ -770,7 +773,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* As of here, particle force variables will be set. */
/* Compute variables required for the force loop */
hydro_prepare_force(p, xp);
hydro_prepare_force(p, xp, cosmo);
/* The particle force values are now set. Do _NOT_
try to read any particle density variables! */
......@@ -1574,6 +1577,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
const int count = c->count;
const int gcount = c->gcount;
const int scount = c->scount;
......@@ -1602,7 +1606,7 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
if (part_is_active(p, e)) {
/* Finish the force loop */
hydro_end_force(p);
hydro_end_force(p, cosmo);
}
}
......
......@@ -34,6 +34,7 @@
/* Local includes. */
#include "active.h"
#include "cell.h"
#include "cosmology.h"
#include "error.h"
#include "gravity.h"
#include "hydro.h"
......@@ -411,7 +412,8 @@ void self_all_force(struct runner *r, struct cell *ci) {
* @brief Compute the force on a single particle brute-force.
*/
void engine_single_density(double *dim, long long int pid,
struct part *restrict parts, int N, int periodic) {
struct part *restrict parts, int N, int periodic,
const struct cosmology *cosmo) {
double r2, dx[3];
float fdx[3];
struct part p;
......@@ -446,7 +448,7 @@ void engine_single_density(double *dim, long long int pid,
}
/* Dump the result. */
hydro_end_density(&p);
hydro_end_density(&p, cosmo);
message("part %lli (h=%e) has wcount=%e, rho=%e.", p.id, p.h,
p.density.wcount, hydro_get_comoving_density(&p));
fflush(stdout);
......
......@@ -48,8 +48,8 @@ velocities\footnote{One additional inconvenience of our choice of
generalised coordinates is that our velocities $\mathbf{v}'$ and
sound-speed $c'$ do not have the same dependencies on the
scale-factor. The signal velocity entering the time-step calculation
will hence read $v_{\rm sig} = a\dot{\mathbf{r}'} + c =
\frac{|\mathbf{v}'|}{a} + a^{1-3(\gamma-1)/2}c'$.}.
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{v}' = a \ddot{r}$. Using the SPH
definition of density, $\rho_i =
\sum_jm_jW(\mathbf{r}_{j}'-\mathbf{r}_{i}',h_i') =
......@@ -109,7 +109,7 @@ viscosity ``tensor'' $\Pi_{ij} \equiv \frac{1}{2}\alpha_{\rm visc} v_{\rm
\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' -
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}'.
......
Supports Markdown
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