diff --git a/src/cosmology.c b/src/cosmology.c
index d304354f3e5068137096c971fed693104924cdb3..b7c1ed2152b063db17a6c2377bac4d98a1587671 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -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);
}
diff --git a/src/cosmology.h b/src/cosmology.h
index 7136b65667195953971060b76ddfd447a5fdf500..4556b039bd0e306dab37a05bc200c3aa2ab8a602 100644
--- a/src/cosmology.h
+++ b/src/cosmology.h
@@ -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);
diff --git a/src/equation_of_state/planetary/hm80.h b/src/equation_of_state/planetary/hm80.h
index 3f9600709836b99c74a4857653214497e0aaa773..5e80c240018756cb57cc8974df4974a6cc53724a 100644
--- a/src/equation_of_state/planetary/hm80.h
+++ b/src/equation_of_state/planetary/hm80.h
@@ -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
diff --git a/src/equation_of_state/planetary/sesame.h b/src/equation_of_state/planetary/sesame.h
index 4deaded44314d6b6204c25a18a6eacc4bdc6ac9a..d958c9b9d09ffe37eefd77ad0384d85bf8c055dd 100644
--- a/src/equation_of_state/planetary/sesame.h
+++ b/src/equation_of_state/planetary/sesame.h
@@ -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;
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 2c3a9c46f0500fb20aa3cfa2e5feb682b3dcec63..237a34283e1b89b06a4de1d6f1890f9a7e1af509 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -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) {}
/**
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 26e3bf97dd1924abbe7380d1eaadce75213344df..4511b2d655b0e7b3293633a466c76757a0237874 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -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 */
diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h
index 7c38896a4bfb0c44a79c509954e11128244e350e..41c870da0cc4630896324b5fdab4b6ca2d4362bc 100644
--- a/src/hydro/GizmoMFM/hydro.h
+++ b/src/hydro/GizmoMFM/hydro.h
@@ -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;
diff --git a/src/hydro/GizmoMFM/hydro_slope_limiters.h b/src/hydro/GizmoMFM/hydro_slope_limiters.h
index 78f2785cdae5dc2334d37e3924dd5b259cca8c05..7c9c759830a4b0ee98412d5a200700c0a148d316 100644
--- a/src/hydro/GizmoMFM/hydro_slope_limiters.h
+++ b/src/hydro/GizmoMFM/hydro_slope_limiters.h
@@ -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
diff --git a/src/hydro/GizmoMFV/hydro.h b/src/hydro/GizmoMFV/hydro.h
index c65f2bb5899b63dbaea4db2f108cf3d914ca78f0..c1a62e3c16b98e98b881f3fe4ddcd539cf842c9d 100644
--- a/src/hydro/GizmoMFV/hydro.h
+++ b/src/hydro/GizmoMFV/hydro.h
@@ -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 */
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 812f8ad72de55ad7990ee6ef88223a401780bc4b..01abe55e7267ca04a7f3e9740c10c681f86f29ea 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -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*/
diff --git a/src/hydro/Planetary/hydro.h b/src/hydro/Planetary/hydro.h
index 134893d9bb01e86de0a17e196d307e6972372c03..106e1a96ae57868c94d1077b74e84909ab0f0830 100644
--- a/src/hydro/Planetary/hydro.h
+++ b/src/hydro/Planetary/hydro.h
@@ -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*/
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index ea086daeeb1e93d7f1476302564fb4182a6fb611..5082db6c792701511972a52fd2fb00a6a45f7271 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -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) {
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index e4b7cf06e083638a94526cc1f9e7212cf19dfad4..3f0f9931ebd557fffcab7e89f3c6297c2fb26474 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -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*/
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index d70d58c6ba508ba4282ac9dd32565478afb40692..cca2a241866fc797055922a48c25cebd6fa1b140 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -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,
diff --git a/src/kick.h b/src/kick.h
index 9d10f1e78d3934c4277c14217cbbc46514e87033..50ecaea498bdd401cc0ac27525ed27986a344c59 100644
--- a/src/kick.h
+++ b/src/kick.h
@@ -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);
}
diff --git a/src/runner.c b/src/runner.c
index 57cfe863d2a1efbe8c62c91d2d4b02d495630c22..6093f1194d24cb6a0af1b2fe31d04be3196930e4 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -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);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
diff --git a/tests/testOutputList.c b/tests/testOutputList.c
index f3aa155577c992caeaee8516a298cb58b94ec2e1..b7df197405ee095cf9bf0a63e8cf7f00585f269f 100644
--- a/tests/testOutputList.c
+++ b/tests/testOutputList.c
@@ -29,12 +29,16 @@
/* Expected values from file */
const double time_values[Ntest] = {
- 0., 10., 12.,
+ 0.,
+ 10.,
+ 12.,
};
/* Expected values from file */
const double a_values[Ntest] = {
- 0.5, 0.1, 0.01,
+ 0.5,
+ 0.1,
+ 0.01,
};
void test_no_cosmo(struct engine *e, char *name, int with_assert) {
@@ -55,7 +59,7 @@ void test_no_cosmo(struct engine *e, char *name, int with_assert) {
output_list_init(&list, e, name, &delta_time, &output_time);
output_list_print(list);
- for(int i = 0; i < Ntest; i++) {
+ for (int i = 0; i < Ntest; i++) {
/* Test last value */
if (with_assert) {
assert(abs(output_time - time_values[i]) < tol);
@@ -69,12 +73,10 @@ void test_no_cosmo(struct engine *e, char *name, int with_assert) {
integertime_t ti_next;
output_list_read_next_time(list, e, name, &ti_next);
- output_time = (double) (ti_next * e->time_base) + e->time_begin;
+ output_time = (double)(ti_next * e->time_base) + e->time_begin;
}
-
output_list_clean(list);
-
};
void test_cosmo(struct engine *e, char *name, int with_assert) {
@@ -93,7 +95,7 @@ void test_cosmo(struct engine *e, char *name, int with_assert) {
output_list_init(&list, e, name, &delta_time, &output_time);
output_list_print(list);
- for(int i = 0; i < Ntest; i++) {
+ for (int i = 0; i < Ntest; i++) {
/* Test last value */
if (with_assert) {
assert(abs(output_time - a_values[i]) < tol);
@@ -107,15 +109,12 @@ void test_cosmo(struct engine *e, char *name, int with_assert) {
integertime_t ti_next;
output_list_read_next_time(list, e, name, &ti_next);
- output_time = (double) exp(ti_next * e->time_base) * e->cosmology->a_begin;
+ output_time = (double)exp(ti_next * e->time_base) * e->cosmology->a_begin;
}
output_list_clean(list);
-
};
-
-
int main(int argc, char *argv[]) {
/* Create a structure to read file into. */
struct swift_params params;
@@ -146,7 +145,7 @@ int main(int argc, char *argv[]) {
int without_assert = 0;
/* Test without cosmo */
test_no_cosmo(&e, "Time", with_assert);
-
+
/* Test with cosmo */
test_cosmo(&e, "Redshift", with_assert);
test_cosmo(&e, "ScaleFactor", with_assert);