diff --git a/src/const.h b/src/const.h
index c290a3e73a524e9003cadb63f8bd45e8b3c51dac..1239d8cd6894aa678f706bc2f67df764e259d533 100644
--- a/src/const.h
+++ b/src/const.h
@@ -55,6 +55,10 @@
#define SLOPE_LIMITER_PER_FACE
#define SLOPE_LIMITER_CELL_WIDE
+/* Options to control the movement of particles for GIZMO_SPH. */
+/* This option disables particle movement */
+//#define GIZMO_FIX_PARTICLES
+
/* Source terms */
#define SOURCETERMS_NONE
//#define SOURCETERMS_SN_FEEDBACK
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index bfb5cd1ce39a9908573c66406f41b56561a870d6..a614d08c30b21f9e7d422bf6b6a09d10d2e89799 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -148,6 +148,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return min(dt_cfl, dt_u_change);
}
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part *p, float dt) {}
+
/**
* @brief Prepares a particle for the density calculation.
*
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 160a2d8b5d25a97cefb2afd5e22d8e6bcea0006e..cc7b422ccbe7c678969df5779a4d4a054c65528e 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -152,6 +152,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return dt_cfl;
}
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part *p, float dt) {}
+
/**
* @brief Prepares a particle for the density calculation.
*
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index c59af05460157a756c15d8ca84af8a7834fde2d3..71a082c6228591a64a4b71b950106670e7513c78 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -1,3 +1,4 @@
+
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
@@ -40,6 +41,27 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return CFL_condition * p->h / fabsf(p->timestepvars.vmax);
}
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * We use this to store the physical time step, since it is used for the flux
+ * exchange during the force loop.
+ *
+ * We also set the active flag of the particle to inactive. It will be set to
+ * active in hydro_init_part, which is called the next time the particle becomes
+ * active.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part* p, float dt) {
+
+ p->force.dt = dt;
+ p->force.active = 0;
+}
+
/**
* @brief Initialises the particles for the first time
*
@@ -58,13 +80,29 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part* p, struct xpart* xp) {
- xp->v_full[0] = p->v[0];
- xp->v_full[1] = p->v[1];
- xp->v_full[2] = p->v[2];
+ const float mass = p->conserved.mass;
p->primitives.v[0] = p->v[0];
p->primitives.v[1] = p->v[1];
p->primitives.v[2] = p->v[2];
+
+ /* we can already initialize the momentum */
+ p->conserved.momentum[0] = mass * p->primitives.v[0];
+ p->conserved.momentum[1] = mass * p->primitives.v[1];
+ p->conserved.momentum[2] = mass * p->primitives.v[2];
+
+ /* and the thermal energy */
+ p->conserved.energy *= mass;
+
+#if defined(GIZMO_FIX_PARTICLES)
+ p->v[0] = 0.;
+ p->v[1] = 0.;
+ p->v[2] = 0.;
+#else
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
+#endif
}
/**
@@ -89,6 +127,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->geometry.matrix_E[2][0] = 0.0f;
p->geometry.matrix_E[2][1] = 0.0f;
p->geometry.matrix_E[2][2] = 0.0f;
+
+ /* Set the active flag to active. */
+ p->force.active = 1;
}
/**
@@ -159,6 +200,12 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->primitives.v[2] = momentum[2] / m;
const float energy = p->conserved.energy;
p->primitives.P = hydro_gamma_minus_one * energy / volume;
+
+ /* sanity checks */
+ if (p->primitives.rho < 0.0f || p->primitives.P < 0.0f) {
+ p->primitives.rho = 0.0f;
+ p->primitives.P = 0.0f;
+ }
}
/**
@@ -180,9 +227,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* restrict p, struct xpart* restrict xp) {
- /* Set the physical time step */
- p->force.dt = get_timestep(p->time_bin, 0.); // MATTHIEU 0
-
/* Initialize time step criterion variables */
p->timestepvars.vmax = 0.0f;
@@ -246,40 +290,14 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* @brief Converts the hydrodynamic variables from the initial condition file to
* conserved variables that can be used during the integration
*
- * Requires the volume to be known.
- *
- * The initial condition file contains a mixture of primitive and conserved
- * variables. Mass is a conserved variable, and we just copy the particle
- * mass into the corresponding conserved quantity. We need the volume to
- * also derive a density, which is then used to convert the internal energy
- * to a pressure. However, we do not actually use these variables anymore.
- * We do need to initialize the linear momentum, based on the mass and the
- * velocity of the particle.
+ * We no longer do this, as the mass needs to be provided in the initial
+ * condition file, and the mass alone is enough to initialize all conserved
+ * variables. This is now done in hydro_first_init_part.
*
* @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part* p, struct xpart* xp) {
-
- const float volume = p->geometry.volume;
- const float m = p->conserved.mass;
- p->primitives.rho = m / volume;
-
- /* first get the initial velocities, as they were overwritten in end_density
- */
- p->primitives.v[0] = p->v[0];
- p->primitives.v[1] = p->v[1];
- p->primitives.v[2] = p->v[2];
-
- p->conserved.momentum[0] = m * p->primitives.v[0];
- p->conserved.momentum[1] = m * p->primitives.v[1];
- p->conserved.momentum[2] = m * p->primitives.v[2];
-
- p->primitives.P =
- hydro_gamma_minus_one * p->conserved.energy * p->primitives.rho;
-
- p->conserved.energy *= m;
-}
+ struct part* p, struct xpart* xp) {}
/**
* @brief Extra operations to be done during the drift
@@ -362,6 +380,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
p->du_dt = 0.0f;
}
+
+#if defined(GIZMO_FIX_PARTICLES)
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+
+ p->du_dt = 0.0f;
+
+ /* disable the smoothing length update, since the smoothing lengths should
+ stay the same for all steps (particles don't move) */
+ p->force.h_dt = 0.0f;
+#endif
}
/**
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index 90448efc7adb8ccecaaa98c7388f89eaa8d16bcd..65f4cb8cf892d1a9948bc043254ddb1581875ed3 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -140,57 +140,61 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
/* time */
- dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
- Wi[2] * pi->primitives.gradients.rho[1] +
- Wi[3] * pi->primitives.gradients.rho[2] +
- Wi[0] * (pi->primitives.gradients.v[0][0] +
- pi->primitives.gradients.v[1][1] +
- pi->primitives.gradients.v[2][2]));
- dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
- Wi[2] * pi->primitives.gradients.v[0][1] +
- Wi[3] * pi->primitives.gradients.v[0][2] +
- pi->primitives.gradients.P[0] / Wi[0]);
- dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
- Wi[2] * pi->primitives.gradients.v[1][1] +
- Wi[3] * pi->primitives.gradients.v[1][2] +
- pi->primitives.gradients.P[1] / Wi[0]);
- dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
- Wi[2] * pi->primitives.gradients.v[2][1] +
- Wi[3] * pi->primitives.gradients.v[2][2] +
- pi->primitives.gradients.P[2] / Wi[0]);
- dWi[4] -=
- 0.5 * mindt * (Wi[1] * pi->primitives.gradients.P[0] +
- Wi[2] * pi->primitives.gradients.P[1] +
- Wi[3] * pi->primitives.gradients.P[2] +
- hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
- pi->primitives.gradients.v[1][1] +
- pi->primitives.gradients.v[2][2]));
-
- dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
- Wj[2] * pj->primitives.gradients.rho[1] +
- Wj[3] * pj->primitives.gradients.rho[2] +
- Wj[0] * (pj->primitives.gradients.v[0][0] +
- pj->primitives.gradients.v[1][1] +
- pj->primitives.gradients.v[2][2]));
- dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
- Wj[2] * pj->primitives.gradients.v[0][1] +
- Wj[3] * pj->primitives.gradients.v[0][2] +
- pj->primitives.gradients.P[0] / Wj[0]);
- dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
- Wj[2] * pj->primitives.gradients.v[1][1] +
- Wj[3] * pj->primitives.gradients.v[1][2] +
- pj->primitives.gradients.P[1] / Wj[0]);
- dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
- Wj[2] * pj->primitives.gradients.v[2][1] +
- Wj[3] * pj->primitives.gradients.v[2][2] +
- pj->primitives.gradients.P[2] / Wj[0]);
- dWj[4] -=
- 0.5 * mindt * (Wj[1] * pj->primitives.gradients.P[0] +
- Wj[2] * pj->primitives.gradients.P[1] +
- Wj[3] * pj->primitives.gradients.P[2] +
- hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
- pj->primitives.gradients.v[1][1] +
- pj->primitives.gradients.v[2][2]));
+ if (Wi[0] > 0.0f) {
+ dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
+ Wi[2] * pi->primitives.gradients.rho[1] +
+ Wi[3] * pi->primitives.gradients.rho[2] +
+ Wi[0] * (pi->primitives.gradients.v[0][0] +
+ pi->primitives.gradients.v[1][1] +
+ pi->primitives.gradients.v[2][2]));
+ dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
+ Wi[2] * pi->primitives.gradients.v[0][1] +
+ Wi[3] * pi->primitives.gradients.v[0][2] +
+ pi->primitives.gradients.P[0] / Wi[0]);
+ dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
+ Wi[2] * pi->primitives.gradients.v[1][1] +
+ Wi[3] * pi->primitives.gradients.v[1][2] +
+ pi->primitives.gradients.P[1] / Wi[0]);
+ dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
+ Wi[2] * pi->primitives.gradients.v[2][1] +
+ Wi[3] * pi->primitives.gradients.v[2][2] +
+ pi->primitives.gradients.P[2] / Wi[0]);
+ dWi[4] -= 0.5 * mindt *
+ (Wi[1] * pi->primitives.gradients.P[0] +
+ Wi[2] * pi->primitives.gradients.P[1] +
+ Wi[3] * pi->primitives.gradients.P[2] +
+ hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
+ pi->primitives.gradients.v[1][1] +
+ pi->primitives.gradients.v[2][2]));
+ }
+
+ if (Wj[0] > 0.0f) {
+ dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
+ Wj[2] * pj->primitives.gradients.rho[1] +
+ Wj[3] * pj->primitives.gradients.rho[2] +
+ Wj[0] * (pj->primitives.gradients.v[0][0] +
+ pj->primitives.gradients.v[1][1] +
+ pj->primitives.gradients.v[2][2]));
+ dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
+ Wj[2] * pj->primitives.gradients.v[0][1] +
+ Wj[3] * pj->primitives.gradients.v[0][2] +
+ pj->primitives.gradients.P[0] / Wj[0]);
+ dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
+ Wj[2] * pj->primitives.gradients.v[1][1] +
+ Wj[3] * pj->primitives.gradients.v[1][2] +
+ pj->primitives.gradients.P[1] / Wj[0]);
+ dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
+ Wj[2] * pj->primitives.gradients.v[2][1] +
+ Wj[3] * pj->primitives.gradients.v[2][2] +
+ pj->primitives.gradients.P[2] / Wj[0]);
+ dWj[4] -= 0.5 * mindt *
+ (Wj[1] * pj->primitives.gradients.P[0] +
+ Wj[2] * pj->primitives.gradients.P[1] +
+ Wj[3] * pj->primitives.gradients.P[2] +
+ hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
+ pj->primitives.gradients.v[1][1] +
+ pj->primitives.gradients.v[2][2]));
+ }
Wi[0] += dWi[0];
Wi[1] += dWi[1];
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index aba6bd53c1c9557929426c11a0986e5f02888874..657119939d95c97ed049df95e9a4133fab135908 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -231,7 +231,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
dtj = pj->force.dt;
/* calculate the maximal signal velocity */
- if (Wi[0] && Wj[0]) {
+ if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
vmax =
sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
} else {
@@ -412,11 +412,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
==> we update particle j if (MODE IS 1) OR (j IS INACTIVE)
*/
- // MATTHIEU
- const integertime_t pj_ti_end = 0; // get_integer_time_end(pj->time_bin);
- const integertime_t pi_ti_end = 0; // get_integer_time_end(pi->time_bin);
-
- if (mode == 1 || pj_ti_end > pi_ti_end) {
+ if (mode == 1 || pj->force.active == 0) {
/* Store mass flux */
mflux = dtj * Anorm * totflux[0];
pj->gravity.mflux[0] -= mflux * dx[0];
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index dcf8cd8fd6cf37c0d7f4ba718746ec7dc3e79b01..818217baee18c2128c8f004e3f2b14ecd454acfc 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -23,6 +23,13 @@
#include "io_properties.h"
#include "riemann.h"
+/* Set the description of the particle movement. */
+#if defined(GIZMO_FIX_PARTICLES)
+#define GIZMO_PARTICLE_MOVEMENT "Fixed particles."
+#else
+#define GIZMO_PARTICLE_MOVEMENT "Particles move with flow velocity."
+#endif
+
/**
* @brief Specifies which particle fields to read from a dataset
*
@@ -155,6 +162,9 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Riemann solver information */
io_write_attribute_s(h_grpsph, "Riemann solver type",
RIEMANN_SOLVER_IMPLEMENTATION);
+
+ /* Particle movement information */
+ io_write_attribute_s(h_grpsph, "Particle movement", GIZMO_PARTICLE_MOVEMENT);
}
/**
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index f6592ca107d8d2c6970f34ebd3929e226b53a355..928011d8201f3f355b00d5fe064267d379278e64 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/Gizmo/hydro_part.h
@@ -178,6 +178,9 @@ struct part {
/* Physical time step of the particle. */
float dt;
+ /* Flag keeping track of whether this is an active or inactive particle. */
+ char active;
+
/* Actual velocity of the particle. */
float v_full[3];
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 20856b7e038855e22aa3776a74ba9f495ff6c93f..56078a82569fb0bc30347d5c01831e9eecfd48f4 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -165,6 +165,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return dt_cfl;
}
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part *p, float dt) {}
+
/**
* @brief Prepares a particle for the density calculation.
*
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index f22bb8a13a8ba4d896a77bd4c4f5e86bed5a5960..20238896f1458d0abebacca4865968a3a671c886 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -152,6 +152,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return dt_cfl;
}
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part *p, float dt) {}
+
/**
* @brief Prepares a particle for the density calculation.
*
diff --git a/src/runner.c b/src/runner.c
index 19c01221d1e261a4a332061350c006f49f3fed79..eb771512ae49ebc51e3ca858bf2bd196d5f787aa 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1096,6 +1096,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
struct xpart *restrict xparts = c->xparts;
struct gpart *restrict gparts = c->gparts;
struct spart *restrict sparts = c->sparts;
+ const double timeBase = e->timeBase;
TIMER_TIC;
@@ -1131,6 +1132,10 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
p->time_bin = get_time_bin(ti_new_step);
if (p->gpart != NULL) p->gpart->time_bin = get_time_bin(ti_new_step);
+ /* Tell the particle what the new physical time step is */
+ float dt = get_timestep(p->time_bin, timeBase);
+ hydro_timestep_extra(p, dt);
+
/* Number of updated particles */
updated++;
if (p->gpart != NULL) g_updated++;
diff --git a/tests/benchmarkInteractions.c b/tests/benchmarkInteractions.c
index 6d6d345bee743d28fb4bdda911bd4bcc4c78205f..e3f558f88dffbab252bf7c06f9e943ff568b6fff 100644
--- a/tests/benchmarkInteractions.c
+++ b/tests/benchmarkInteractions.c
@@ -91,7 +91,10 @@ struct part *make_particles(size_t count, double *offset, double spacing,
p->h = h;
p->id = ++(*partId);
+
+#if !defined(GIZMO_SPH)
p->mass = 1.0f;
+#endif
/* Place rest of particles around the test particle
* with random position within a unit sphere. */
@@ -110,7 +113,9 @@ struct part *make_particles(size_t count, double *offset, double spacing,
p->h = h;
p->id = ++(*partId);
+#if !defined(GIZMO_SPH)
p->mass = 1.0f;
+#endif
}
return particles;
}
@@ -120,6 +125,7 @@ struct part *make_particles(size_t count, double *offset, double spacing,
*/
void prepare_force(struct part *parts, size_t count) {
+#if !defined(GIZMO_SPH)
struct part *p;
for (size_t i = 0; i < count; ++i) {
p = &parts[i];
@@ -130,6 +136,7 @@ void prepare_force(struct part *parts, size_t count) {
p->force.v_sig = 0.0f;
p->force.h_dt = 0.0f;
}
+#endif
}
/**
@@ -144,10 +151,19 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
"%13e %13e %13e %13e "
"%13e %13e %13e\n",
p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2],
- p->a_hydro[0], p->a_hydro[1], p->a_hydro[2], p->rho,
- p->density.rho_dh, p->density.wcount, p->density.wcount_dh,
- p->force.h_dt, p->force.v_sig,
-#if defined(MINIMAL_SPH)
+ p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
+#if defined(GIZMO_SPH)
+ 0., 0.,
+#else
+ p->rho, p->density.rho_dh,
+#endif
+ p->density.wcount, p->density.wcount_dh, p->force.h_dt,
+#if defined(GIZMO_SPH)
+ 0.,
+#else
+ p->force.v_sig,
+#endif
+#if defined(MINIMAL_SPH) || defined(GIZMO_SPH)
0., 0., 0., 0.
#else
p->density.div_v, p->density.rot_v[0], p->density.rot_v[1],