Commit 95ef3dbe authored by Bert Vandenbroucke's avatar Bert Vandenbroucke Committed by Matthieu Schaller
Browse files

Gizmo cosmology

parent c49a9670
......@@ -265,6 +265,29 @@ double hydro_kick_integrand(double a, void *param) {
return (1. / H) * pow(a_inv, 3. * hydro_gamma_minus_one) * a_inv;
}
/**
* @brief Computes \f$a dt\f$ for the current cosmology.
*
* @param a The scale-factor of interest.
* @param param The current #cosmology.
*/
double hydro_kick_corr_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 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;
}
/**
* @brief Computes \f$ dt \f$ for the current cosmology.
*
......@@ -306,6 +329,8 @@ void cosmology_init_tables(struct cosmology *c) {
(double *)malloc(cosmology_table_length * sizeof(double));
c->hydro_kick_fac_interp_table =
(double *)malloc(cosmology_table_length * sizeof(double));
c->hydro_kick_corr_interp_table =
(double *)malloc(cosmology_table_length * sizeof(double));
c->time_interp_table =
(double *)malloc(cosmology_table_length * sizeof(double));
c->scale_factor_interp_table =
......@@ -354,6 +379,16 @@ void cosmology_init_tables(struct cosmology *c) {
c->hydro_kick_fac_interp_table[i] = result;
}
/* Integrate the kick correction factor \int_{a_begin}^{a_table[i]} a dt */
F.function = &hydro_kick_corr_integrand;
for (int i = 0; i < cosmology_table_length; i++) {
gsl_integration_qag(&F, a_begin, a_table[i], 0, 1.0e-10, GSL_workspace_size,
GSL_INTEG_GAUSS61, space, &result, &abserr);
/* Store result */
c->hydro_kick_corr_interp_table[i] = result;
}
/* Integrate the time \int_{a_begin}^{a_table[i]} dt */
F.function = &time_integrand;
for (int i = 0; i < cosmology_table_length; i++) {
......@@ -374,7 +409,6 @@ void cosmology_init_tables(struct cosmology *c) {
GSL_INTEG_GAUSS61, space, &result, &abserr);
c->universe_age_at_present_day = result;
/* Update the times */
c->time_begin = cosmology_get_time_since_big_bang(c, c->a_begin);
c->time_end = cosmology_get_time_since_big_bang(c, c->a_end);
......@@ -383,8 +417,7 @@ void cosmology_init_tables(struct cosmology *c) {
* Inverse t(a)
*/
const double delta_t =
(c->time_end - c->time_begin) / cosmology_table_length;
const double delta_t = (c->time_end - c->time_begin) / cosmology_table_length;
/* index in the time_interp_table */
int i_a = 0;
......@@ -411,7 +444,7 @@ void cosmology_init_tables(struct cosmology *c) {
/* Compute interpolated scale factor */
double log_a = c->log_a_begin + scale * (c->log_a_end - c->log_a_begin) /
cosmology_table_length;
cosmology_table_length;
c->scale_factor_interp_table[i_time] = exp(log_a) - c->a_begin;
}
......@@ -625,6 +658,35 @@ double cosmology_get_hydro_kick_factor(const struct cosmology *c,
return int_end - int_start;
}
/**
* @brief Computes the cosmology factor that enters the hydro kick correction
* operator for the meshless schemes (GIZMO-MFV).
*
* Computes \f$ \int_{a_start}^{a_end} a dt \f$ using the interpolation table.
*
* @param c The current #cosmology.
* @param ti_start the (integer) time of the start of the drift.
* @param ti_end the (integer) time of the end of the drift.
*/
double cosmology_get_corr_kick_factor(const struct cosmology *c,
integertime_t ti_start,
integertime_t ti_end) {
#ifdef SWIFT_DEBUG_CHECKS
if (ti_end < ti_start) error("ti_end must be >= ti_start");
#endif
const double a_start = c->log_a_begin + ti_start * c->time_base;
const double a_end = c->log_a_begin + ti_end * c->time_base;
const double int_start = interp_table(c->hydro_kick_corr_interp_table,
a_start, c->log_a_begin, c->log_a_end);
const double int_end = interp_table(c->hydro_kick_corr_interp_table, a_end,
c->log_a_begin, c->log_a_end);
return int_end - int_start;
}
/**
* @brief Computes the cosmology factor that enters the thermal variable kick
* operator.
......@@ -718,6 +780,7 @@ void cosmology_clean(struct cosmology *c) {
free(c->drift_fac_interp_table);
free(c->grav_kick_fac_interp_table);
free(c->hydro_kick_fac_interp_table);
free(c->hydro_kick_corr_interp_table);
free(c->time_interp_table);
free(c->scale_factor_interp_table);
}
......
......@@ -153,6 +153,9 @@ struct cosmology {
/*! Kick factor (hydro) interpolation table */
double *hydro_kick_fac_interp_table;
/*! Kick factor (hydro correction) interpolation table (GIZMO-MFV only) */
double *hydro_kick_corr_interp_table;
/*! Time interpolation table */
double *time_interp_table;
......@@ -180,6 +183,9 @@ double cosmology_get_hydro_kick_factor(const struct cosmology *cosmo,
double cosmology_get_therm_kick_factor(const struct cosmology *cosmo,
integertime_t ti_start,
integertime_t ti_end);
double cosmology_get_corr_kick_factor(const struct cosmology *cosmo,
integertime_t ti_start,
integertime_t ti_end);
double cosmology_get_delta_time(const struct cosmology *c,
integertime_t ti_start, integertime_t ti_end);
......
......@@ -91,7 +91,8 @@ INLINE static void load_table_HM80(struct HM80_params *mat, char *table_file) {
// Ignore header lines
char buffer[100];
for (int i = 0; i < 4; i++) {
fgets(buffer, 100, f);
if (fgets(buffer, 100, f) == NULL)
error("Something incorrect happening with the file header.");
}
// Table properties
......
......@@ -87,7 +87,8 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat,
// Ignore header lines
char buffer[100];
for (int i = 0; i < 5; i++) {
fgets(buffer, 100, f);
if (fgets(buffer, 100, f) == NULL)
error("Something incorrect happening with the file header.");
}
float ignore;
......
......@@ -506,6 +506,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt_therm,
float dt_grav, float dt_hydro, float dt_kick_corr,
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {}
/**
......
......@@ -530,6 +530,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt_therm,
float dt_grav, float dt_hydro, float dt_kick_corr,
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* Do not decrease the entropy by more than a factor of 2 */
......
......@@ -60,10 +60,9 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
sqrtf(hydro_gamma * p->P / p->rho);
vmax = max(vmax, p->timestepvars.vmax);
// MATTHIEU: Bert is this correct? Do we need more cosmology terms here?
const float psize =
cosmo->a * powf(p->geometry.volume / hydro_dimension_unit_sphere,
hydro_dimension_inv);
const float psize = cosmo->a * cosmo->a *
powf(p->geometry.volume / hydro_dimension_unit_sphere,
hydro_dimension_inv);
float dt = FLT_MAX;
if (vmax > 0.0f) {
dt = psize / vmax;
......@@ -502,7 +501,10 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part* p, struct xpart* xp, const struct cosmology* cosmo) {}
struct part* p, struct xpart* xp, const struct cosmology* cosmo) {
p->conserved.energy /= cosmo->a_factor_internal_energy;
}
/**
* @brief Extra operations to be done during the drift
......@@ -541,23 +543,17 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
#if !defined(EOS_ISOTHERMAL_GAS)
#ifdef GIZMO_TOTAL_ENERGY
const float Etot = p->conserved.energy + p->flux.energy * dt_therm;
const float Etot = p->conserved.energy + p->flux.energy * dt_drift;
const float v2 =
(p->v[0] * p->v[0] + p->v[1] * p->v[1] + p->v[2] * p->v[2]);
const float u = (Etot * m_inv - 0.5f * v2);
#else
const float u = (p->conserved.energy + p->flux.energy * dt_therm) * m_inv;
const float u = (p->conserved.energy + p->flux.energy * dt_drift) * m_inv;
#endif
p->P = hydro_gamma_minus_one * u * p->rho;
#endif
}
/* we use a sneaky way to get the gravitational contribtuion to the
velocity update */
p->v[0] += p->v[0] - xp->v_full[0];
p->v[1] += p->v[1] - xp->v_full[1];
p->v[2] += p->v[2] - xp->v_full[2];
#ifdef SWIFT_DEBUG_CHECKS
if (p->h <= 0.0f) {
error("Zero or negative smoothing length (%g)!", p->h);
......@@ -595,29 +591,40 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
/**
* @brief Extra operations done during the kick
*
* Not used for GIZMO.
*
* @param p Particle to act upon.
* @param xp Extended particle data to act upon.
* @param dt Physical time step.
* @param half_dt Half the physical time step.
* @param dt_therm Thermal energy time-step @f$\frac{dt}{a^2}@f$.
* @param dt_grav Gravity time-step @f$\frac{dt}{a}@f$.
* @param dt_hydro Hydro acceleration time-step
* @f$\frac{dt}{a^{3(\gamma{}-1)}}@f$.
* @param dt_kick_corr Gravity correction time-step @f$adt@f$.
* @param cosmo Cosmology.
* @param hydro_props Additional hydro properties.
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt, const struct cosmology* cosmo,
struct part* p, struct xpart* xp, float dt_therm, float dt_grav,
float dt_hydro, float dt_kick_corr, const struct cosmology* cosmo,
const struct hydro_props* hydro_props) {
float a_grav[3];
/* Update conserved variables (note: the mass does not change). */
p->conserved.momentum[0] += p->flux.momentum[0] * dt;
p->conserved.momentum[1] += p->flux.momentum[1] * dt;
p->conserved.momentum[2] += p->flux.momentum[2] * dt;
p->conserved.momentum[0] += p->flux.momentum[0] * dt_therm;
p->conserved.momentum[1] += p->flux.momentum[1] * dt_therm;
p->conserved.momentum[2] += p->flux.momentum[2] * dt_therm;
#if defined(EOS_ISOTHERMAL_GAS)
/* We use the EoS equation in a sneaky way here just to get the constant u */
p->conserved.energy =
p->conserved.mass * gas_internal_energy_from_entropy(0.0f, 0.0f);
#else
p->conserved.energy += p->flux.energy * dt;
p->conserved.energy += p->flux.energy * dt_therm;
#endif
#ifndef HYDRO_GAMMA_5_3
const float Pcorr = (dt_hydro - dt_therm) * p->geometry.volume;
p->conserved.momentum[0] -= Pcorr * p->gradients.P[0];
p->conserved.momentum[1] -= Pcorr * p->gradients.P[1];
p->conserved.momentum[2] -= Pcorr * p->gradients.P[2];
#endif
/* Apply the minimal energy limit */
......@@ -655,9 +662,9 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* Kick the momentum for half a time step */
/* Note that this also affects the particle movement, as the velocity for
the particles is set after this. */
p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0];
p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
p->conserved.momentum[0] += dt_grav * p->conserved.mass * a_grav[0];
p->conserved.momentum[1] += dt_grav * p->conserved.mass * a_grav[1];
p->conserved.momentum[2] += dt_grav * p->conserved.mass * a_grav[2];
}
/* Set the velocities: */
......@@ -844,7 +851,6 @@ __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
v[2] = p->v[2];
}
// MATTHIEU: Bert is this correct?
v[0] += xp->a_grav[0] * dt_kick_grav;
v[1] += xp->a_grav[1] * dt_kick_grav;
v[2] += xp->a_grav[2] * dt_kick_grav;
......
......@@ -47,8 +47,8 @@
* @param r Distance between particle i and particle j.
*/
__attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
float *Wi, float *Wj, float *dWi, float *dWj, float *xij_i, float *xij_j,
float r) {}
float *Wi, float *Wj, float *dWi, float *dWj, const float *xij_i,
const float *xij_j, float r) {}
#endif
......
......@@ -68,14 +68,13 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
vmax = max(vmax, p->timestepvars.vmax);
// MATTHIEU: Bert is this correct? Do we need more cosmology terms here?
const float psize =
cosmo->a * powf(p->geometry.volume / hydro_dimension_unit_sphere,
hydro_dimension_inv);
const float psize = powf(p->geometry.volume / hydro_dimension_unit_sphere,
hydro_dimension_inv);
float dt = FLT_MAX;
if (vmax > 0.) {
dt = psize / vmax;
}
return CFL_condition * dt;
return cosmo->a * cosmo->a * CFL_condition * dt;
}
/**
......@@ -544,7 +543,10 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part* p, struct xpart* xp, const struct cosmology* cosmo) {}
struct part* p, struct xpart* xp, const struct cosmology* cosmo) {
p->conserved.energy /= cosmo->a_factor_internal_energy;
}
/**
* @brief Extra operations to be done during the drift
......@@ -579,16 +581,16 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* drift the primitive variables based on the old fluxes */
if (p->geometry.volume > 0.) {
p->primitives.rho += p->conserved.flux.mass * dt_drift / p->geometry.volume;
p->primitives.rho += p->conserved.flux.mass * dt_therm / p->geometry.volume;
}
if (p->conserved.mass > 0.) {
p->primitives.v[0] +=
p->conserved.flux.momentum[0] * dt_drift / p->conserved.mass;
p->conserved.flux.momentum[0] * dt_therm / p->conserved.mass;
p->primitives.v[1] +=
p->conserved.flux.momentum[1] * dt_drift / p->conserved.mass;
p->conserved.flux.momentum[1] * dt_therm / p->conserved.mass;
p->primitives.v[2] +=
p->conserved.flux.momentum[2] * dt_drift / p->conserved.mass;
p->conserved.flux.momentum[2] * dt_therm / p->conserved.mass;
#if !defined(EOS_ISOTHERMAL_GAS)
#ifdef GIZMO_TOTAL_ENERGY
......@@ -607,7 +609,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
#endif
}
/* we use a sneaky way to get the gravitational contribtuion to the
/* we use a sneaky way to get the gravitational contribution to the
velocity update */
p->primitives.v[0] += p->v[0] - xp->v_full[0];
p->primitives.v[1] += p->v[1] - xp->v_full[1];
......@@ -649,15 +651,19 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
/**
* @brief Extra operations done during the kick
*
* Not used for GIZMO.
*
* @param p Particle to act upon.
* @param xp Extended particle data to act upon.
* @param dt Physical time step.
* @param half_dt Half the physical time step.
* @param dt_therm Thermal energy time-step @f$\frac{dt}{a^2}@f$.
* @param dt_grav Gravity time-step @f$\frac{dt}{a}@f$.
* @param dt_hydro Hydro acceleration time-step
* @f$\frac{dt}{a^{3(\gamma{}-1)}}@f$.
* @param dt_kick_corr Gravity correction time-step @f$adt@f$.
* @param cosmo Cosmology.
* @param hydro_props Additional hydro properties.
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt, const struct cosmology* cosmo,
struct part* p, struct xpart* xp, float dt_therm, float dt_grav,
float dt_hydro, float dt_kick_corr, const struct cosmology* cosmo,
const struct hydro_props* hydro_props) {
float a_grav[3];
......@@ -672,35 +678,48 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
a_grav[2] = p->gpart->a_grav[2];
#ifdef GIZMO_TOTAL_ENERGY
p->conserved.energy += dt * (p->conserved.momentum[0] * a_grav[0] +
p->conserved.momentum[1] * a_grav[1] +
p->conserved.momentum[2] * a_grav[2]);
p->conserved.energy += dt_grav * (p->conserved.momentum[0] * a_grav[0] +
p->conserved.momentum[1] * a_grav[1] +
p->conserved.momentum[2] * a_grav[2]);
#endif
/* Kick the momentum for half a time step */
/* Note that this also affects the particle movement, as the velocity for
the particles is set after this. */
p->conserved.momentum[0] += p->conserved.mass * a_grav[0] * dt;
p->conserved.momentum[1] += p->conserved.mass * a_grav[1] * dt;
p->conserved.momentum[2] += p->conserved.mass * a_grav[2] * dt;
p->conserved.momentum[0] += p->conserved.mass * a_grav[0] * dt_grav;
p->conserved.momentum[1] += p->conserved.mass * a_grav[1] * dt_grav;
p->conserved.momentum[2] += p->conserved.mass * a_grav[2] * dt_grav;
p->conserved.energy -=
0.5f * dt *
0.5f * dt_kick_corr *
(p->gravity.mflux[0] * a_grav[0] + p->gravity.mflux[1] * a_grav[1] +
p->gravity.mflux[2] * a_grav[2]);
}
/* Update conserved variables. */
p->conserved.mass += p->conserved.flux.mass * dt;
p->conserved.momentum[0] += p->conserved.flux.momentum[0] * dt;
p->conserved.momentum[1] += p->conserved.flux.momentum[1] * dt;
p->conserved.momentum[2] += p->conserved.flux.momentum[2] * dt;
p->conserved.mass += p->conserved.flux.mass * dt_therm;
p->conserved.momentum[0] += p->conserved.flux.momentum[0] * dt_therm;
p->conserved.momentum[1] += p->conserved.flux.momentum[1] * dt_therm;
p->conserved.momentum[2] += p->conserved.flux.momentum[2] * dt_therm;
#if defined(EOS_ISOTHERMAL_GAS)
/* We use the EoS equation in a sneaky way here just to get the constant u */
p->conserved.energy =
p->conserved.mass * gas_internal_energy_from_entropy(0.f, 0.f);
#else
p->conserved.energy += p->conserved.flux.energy * dt;
p->conserved.energy += p->conserved.flux.energy * dt_therm;
#endif
#ifndef HYDRO_GAMMA_5_3
const float Pcorr = (dt_hydro - dt_therm) * p->geometry.volume;
p->conserved.momentum[0] -= Pcorr * p->primitives.gradients.P[0];
p->conserved.momentum[1] -= Pcorr * p->primitives.gradients.P[1];
p->conserved.momentum[2] -= Pcorr * p->primitives.gradients.P[2];
#ifdef GIZMO_TOTAL_ENERGY
p->conserved.energy -=
Pcorr * (p->primitives.v[0] * p->primitives.gradients.P[0] +
p->primitives.v[1] * p->primitives.gradients.P[1] +
p->primitives.v[2] * p->primitives.gradients.P[2]);
#endif
#endif
/* Apply the minimal energy limit */
......
......@@ -527,6 +527,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt_therm,
float dt_grav, float dt_hydro, float dt_kick_corr,
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* Do not decrease the energy by more than a factor of 2*/
......
......@@ -568,6 +568,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt_therm,
float dt_grav, float dt_hydro, float dt_kick_corr,
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* Do not decrease the energy by more than a factor of 2*/
......
......@@ -581,6 +581,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt_therm,
float dt_grav, float dt_hydro, float dt_kick_corr,
const struct cosmology *cosmo,
const struct hydro_props *restrict hydro_properties) {
......
......@@ -541,6 +541,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt_therm,
float dt_grav, float dt_hydro, float dt_kick_corr,
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* Do not decrease the entropy (temperature) by more than a factor of 2*/
......
......@@ -443,7 +443,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param dt Physical time step.
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt, const struct cosmology* cosmo,
struct part* p, struct xpart* xp, float dt, float dt_grav, float dt_hydro,
float dt_kick_corr, const struct cosmology* cosmo,
const struct hydro_props* hydro_props) {
/* Update the conserved variables. We do this here and not in the kick,
......
......@@ -68,6 +68,7 @@ __attribute__((always_inline)) INLINE static void kick_gpart(
* @param dt_kick_hydro The kick time-step for hydro accelerations.
* @param dt_kick_grav The kick time-step for gravity accelerations.
* @param dt_kick_therm The kick time-step for changes in thermal state.
* @param dt_kick_corr The kick time-step for the gizmo-mfv gravity correction.
* @param cosmo The cosmological model.
* @param hydro_props The constants used in the scheme
* @param ti_start The starting (integer) time of the kick (for debugging
......@@ -76,9 +77,9 @@ __attribute__((always_inline)) INLINE static void kick_gpart(
*/
__attribute__((always_inline)) INLINE static void kick_part(
struct part *restrict p, struct xpart *restrict xp, double dt_kick_hydro,
double dt_kick_grav, double dt_kick_therm, const struct cosmology *cosmo,
const struct hydro_props *hydro_props, integertime_t ti_start,
integertime_t ti_end) {
double dt_kick_grav, double dt_kick_therm, double dt_kick_corr,
const struct cosmology *cosmo, const struct hydro_props *hydro_props,
integertime_t ti_start, integertime_t ti_end) {
#ifdef SWIFT_DEBUG_CHECKS
if (p->ti_kick != ti_start)
......@@ -110,7 +111,8 @@ __attribute__((always_inline)) INLINE static void kick_part(
}
/* Extra kick work */
hydro_kick_extra(p, xp, dt_kick_therm, cosmo, hydro_props);
hydro_kick_extra(p, xp, dt_kick_therm, dt_kick_grav, dt_kick_hydro,
dt_kick_corr, cosmo, hydro_props);
if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt_kick_grav);
}
......
......@@ -1142,7 +1142,7 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
#endif
/* Time interval for this half-kick */
double dt_kick_grav, dt_kick_hydro, dt_kick_therm;
double dt_kick_grav, dt_kick_hydro, dt_kick_therm, dt_kick_corr;
if (with_cosmology) {
dt_kick_hydro = cosmology_get_hydro_kick_factor(
cosmo, ti_begin, ti_begin + ti_step / 2);
......@@ -1150,15 +1150,19 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
ti_begin + ti_step / 2);
dt_kick_therm = cosmology_get_therm_kick_factor(
cosmo, ti_begin, ti_begin + ti_step / 2);
dt_kick_corr = cosmology_get_corr_kick_factor(cosmo, ti_begin,
ti_begin + ti_step / 2);
} else {
dt_kick_hydro = (ti_step / 2) * time_base;
dt_kick_grav = (ti_step / 2) * time_base;
dt_kick_therm = (ti_step / 2) * time_base;
dt_kick_corr = (ti_step / 2) * time_base;
}
/* do the kick */
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, cosmo,
hydro_props, ti_begin, ti_begin + ti_step / 2);
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm,
dt_kick_corr, cosmo, hydro_props, ti_begin,
ti_begin + ti_step / 2);
/* Update the accelerations to be used in the drift for hydro */
if (p->gpart != NULL) {
......@@ -1308,7 +1312,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
ti_begin, ti_step, p->time_bin, ti_current);
#endif
/* Time interval for this half-kick */
double dt_kick_grav, dt_kick_hydro, dt_kick_therm;
double dt_kick_grav, dt_kick_hydro, dt_kick_therm, dt_kick_corr;
if (with_cosmology) {
dt_kick_hydro = cosmology_get_hydro_kick_factor(
cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
......@@ -1316,15 +1320,19 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
dt_kick_therm = cosmology_get_therm_kick_factor(
cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
dt_kick_corr = cosmology_get_corr_kick_factor(
cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
} else {
dt_kick_hydro = (ti_step / 2) * time_base;
dt_kick_grav = (ti_step / 2) * time_base;
dt_kick_therm = (ti_step / 2) * time_base;
dt_kick_corr = (ti_step / 2) * time_base;
}
/* Finish the time-step with a second half-kick */
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, cosmo,
hydro_props, ti_begin + ti_step / 2, ti_begin + ti_step);
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm,
dt_kick_corr, cosmo, hydro_props, ti_begin + ti_step / 2,
ti_begin + ti_step);