diff --git a/examples/SodShock_1D/sodShock.yml b/examples/SodShock_1D/sodShock.yml
index a5759109e48b7e7eb9cbd15957cf438edd909f1f..a3108f4adcb3c19875c5229e2edebb915a5b1e19 100644
--- a/examples/SodShock_1D/sodShock.yml
+++ b/examples/SodShock_1D/sodShock.yml
@@ -9,15 +9,15 @@ InternalUnitSystem:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 0.2 # The end time of the simulation (in internal units).
+ time_end: 0.1 # The end time of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+ dt_max: 0.01 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: sodShock # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
- delta_time: 0.2 # Time difference between consecutive outputs (in internal units)
+ delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/main.c b/examples/main.c
index 7f3829a7a7ed6379c8b307d211ecbe9c2911d7b6..9ab40413b048cf08e035aeec6775bdc8ab1ab9ca 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -386,7 +386,7 @@ int main(int argc, char *argv[]) {
struct gravity_props gravity_properties;
if (with_self_gravity) gravity_props_init(&gravity_properties, params);
-#if defined(SHADOWSWIFT)
+#if defined(SHADOWFAX_SPH)
/* Override the variables governing the density iteration
(see src/hydro/Shadowswift/hydro.h for full explanation) */
hydro_properties.target_neighbours = 1.0f;
diff --git a/src/const.h b/src/const.h
index 88c1a1af1cfc36cc401fdfea0b077f79fcd13bc0..322421c8da3bf1811a1b15fd483dd5a524c9f7fe 100644
--- a/src/const.h
+++ b/src/const.h
@@ -57,6 +57,10 @@
//#define GIZMO_FIX_PARTICLES
//#define GIZMO_TOTAL_ENERGY
+/* Options to control SHADOWFAX_SPH */
+/* This option disables cell movement */
+//#define SHADOWFAX_FIX_CELLS
+
/* Source terms */
#define SOURCETERMS_NONE
//#define SOURCETERMS_SN_FEEDBACK
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 25d8b30f8b643c7357286f5df7ea1414648160b7..1b0430a31ef6b1d13c4445dae5f92a6444bb6d29 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -57,7 +57,11 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* @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) {}
+ struct part* p, float dt) {
+
+ p->force.dt = dt;
+ p->force.active = 0;
+}
/**
* @brief Initialises the particles for the first time
@@ -77,13 +81,27 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
__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];
+
+ 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];
+
+ p->conserved.energy *= mass;
+
+#if defined(SHADOWFAX_FIX_CELLS)
+ p->v[0] = 0.;
+ p->v[1] = 0.;
+ p->v[2] = 0.;
+#endif
+
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
}
/**
@@ -102,6 +120,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->rho = 0.0f;
voronoi_cell_init(&p->cell, p->x);
+
+ /* Set the active flag to active. */
+ p->force.active = 1;
}
/**
@@ -256,25 +277,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* @param xp The extended particle data to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part* p, struct xpart* xp) {
-
- float volume;
- float m;
- float momentum[3];
- volume = p->cell.volume;
-
- p->mass = m = p->conserved.mass;
- p->primitives.rho = m / volume;
-
- p->conserved.momentum[0] = momentum[0] = m * p->primitives.v[0];
- p->conserved.momentum[1] = momentum[1] = m * p->primitives.v[1];
- p->conserved.momentum[2] = 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
@@ -302,8 +305,19 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
* @param p Particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_end_force(
- struct part* p) {
+ struct part* p) {}
+/**
+ * @brief Extra operations done during the kick
+ *
+ * Not used for Shadowswift.
+ *
+ * @param p Particle to act upon.
+ * @param xp Extended particle data to act upon.
+ * @param dt Physical time step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+ struct part* p, struct xpart* xp, float dt) {
float vcell[3];
/* Update the conserved variables. We do this here and not in the kick,
@@ -323,69 +337,58 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
p->conserved.flux.momentum[2] = 0.0f;
p->conserved.flux.energy = 0.0f;
- /* Set the hydro acceleration, based on the new momentum and mass */
- /* NOTE: the momentum and mass are only correct for active particles, since
- only active particles have received flux contributions from all their
- neighbours. Since this method is only called for active particles,
- this is indeed the case. */
- if (p->force.dt) {
- /* We want the cell velocity to be as close as possible to the fluid
- velocity */
- vcell[0] = p->conserved.momentum[0] / p->conserved.mass;
- vcell[1] = p->conserved.momentum[1] / p->conserved.mass;
- vcell[2] = p->conserved.momentum[2] / p->conserved.mass;
-
- /* To prevent stupid things like cell crossovers or generators that move
- outside their cell, we steer the motion of the cell somewhat */
- if (p->primitives.rho) {
- float centroid[3], d[3];
- float volume, csnd, R, vfac, fac, dnrm;
- voronoi_get_centroid(&p->cell, centroid);
- d[0] = centroid[0] - p->x[0];
- d[1] = centroid[1] - p->x[1];
- d[2] = centroid[2] - p->x[2];
- dnrm = sqrtf(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);
- csnd = sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
- volume = p->cell.volume;
- R = get_radius_dimension_sphere(volume);
- fac = 4.0f * dnrm / R;
- if (fac > 0.9f) {
- if (fac < 1.1f) {
- vfac = csnd * (dnrm - 0.225f * R) / dnrm / (0.05f * R);
- } else {
- vfac = csnd / dnrm;
- }
+ /* We want the cell velocity to be as close as possible to the fluid
+ velocity */
+ vcell[0] = p->conserved.momentum[0] / p->conserved.mass;
+ vcell[1] = p->conserved.momentum[1] / p->conserved.mass;
+ vcell[2] = p->conserved.momentum[2] / p->conserved.mass;
+
+ /* To prevent stupid things like cell crossovers or generators that move
+ outside their cell, we steer the motion of the cell somewhat */
+ if (p->primitives.rho) {
+ float centroid[3], d[3];
+ float volume, csnd, R, vfac, fac, dnrm;
+ voronoi_get_centroid(&p->cell, centroid);
+ d[0] = centroid[0] - p->x[0];
+ d[1] = centroid[1] - p->x[1];
+ d[2] = centroid[2] - p->x[2];
+ dnrm = sqrtf(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);
+ csnd = sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+ volume = p->cell.volume;
+ R = get_radius_dimension_sphere(volume);
+ fac = 4.0f * dnrm / R;
+ if (fac > 0.9f) {
+ if (fac < 1.1f) {
+ vfac = csnd * (dnrm - 0.225f * R) / dnrm / (0.05f * R);
} else {
- vfac = 0.0f;
+ vfac = csnd / dnrm;
}
- vcell[0] += vfac * d[0];
- vcell[1] += vfac * d[1];
- vcell[2] += vfac * d[2];
+ } else {
+ vfac = 0.0f;
}
-
- /* We know the desired velocity; now make sure the particle is accelerated
- accordingly during the next time step */
- p->a_hydro[0] = (vcell[0] - p->force.v_full[0]) / p->force.dt;
- p->a_hydro[1] = (vcell[1] - p->force.v_full[1]) / p->force.dt;
- p->a_hydro[2] = (vcell[2] - p->force.v_full[2]) / p->force.dt;
- } else {
- p->a_hydro[0] = 0.0f;
- p->a_hydro[1] = 0.0f;
- p->a_hydro[2] = 0.0f;
+ vcell[0] += vfac * d[0];
+ vcell[1] += vfac * d[1];
+ vcell[2] += vfac * d[2];
}
-}
-/**
- * @brief Extra operations done during the kick
- *
- * Not used for Shadowswift.
- *
- * @param p Particle to act upon.
- * @param xp Extended particle data to act upon.
- * @param dt Physical time step.
- */
-__attribute__((always_inline)) INLINE static void hydro_kick_extra(
- struct part* p, struct xpart* xp, float dt) {}
+#if defined(SHADOWFAX_FIX_CELLS)
+ xp->v_full[0] = 0.;
+ xp->v_full[1] = 0.;
+ xp->v_full[2] = 0.;
+
+ p->v[0] = 0.;
+ p->v[1] = 0.;
+ p->v[2] = 0.;
+#else
+ xp->v_full[0] = vcell[0];
+ xp->v_full[1] = vcell[1];
+ xp->v_full[2] = vcell[2];
+
+ p->v[0] = xp->v_full[0];
+ p->v[1] = xp->v_full[1];
+ p->v[2] = xp->v_full[2];
+#endif
+}
/**
* @brief Returns the internal energy of a particle
diff --git a/src/hydro/Shadowswift/hydro_debug.h b/src/hydro/Shadowswift/hydro_debug.h
index bc43fda1e1a1caba6b47a661beaf1c030f6089cc..cd62f43b9e5d54a8241699c183e8d59c0d77d27f 100644
--- a/src/hydro/Shadowswift/hydro_debug.h
+++ b/src/hydro/Shadowswift/hydro_debug.h
@@ -24,8 +24,6 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"v=[%.3e,%.3e,%.3e], "
"a=[%.3e,%.3e,%.3e], "
"h=%.3e, "
- "ti_begin=%d, "
- "ti_end=%d, "
"primitives={"
"v=[%.3e,%.3e,%.3e], "
"rho=%.3e, "
@@ -50,9 +48,9 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"wcount=%.3e}, "
"mass=%.3e\n",
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->h, p->ti_begin, p->ti_end,
- p->primitives.v[0], p->primitives.v[1], p->primitives.v[2],
- p->primitives.rho, p->primitives.P, p->primitives.gradients.rho[0],
+ p->a_hydro[1], p->a_hydro[2], p->h, p->primitives.v[0],
+ p->primitives.v[1], p->primitives.v[2], p->primitives.rho,
+ p->primitives.P, p->primitives.gradients.rho[0],
p->primitives.gradients.rho[1], p->primitives.gradients.rho[2],
p->primitives.gradients.v[0][0], p->primitives.gradients.v[0][1],
p->primitives.gradients.v[0][2], p->primitives.gradients.v[1][0],
diff --git a/src/hydro/Shadowswift/hydro_iact.h b/src/hydro/Shadowswift/hydro_iact.h
index 8997155a0cb30848a74b220811970421b6d3396a..f99f8f029fa634206f03ed9854acf395264848bb 100644
--- a/src/hydro/Shadowswift/hydro_iact.h
+++ b/src/hydro/Shadowswift/hydro_iact.h
@@ -286,7 +286,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
UPDATE particle j.
==> we update particle j if (MODE IS 1) OR (j IS INACTIVE)
*/
- if (mode == 1 || pj->ti_end > pi->ti_end) {
+ if (mode == 1 || pj->force.active == 0) {
pj->conserved.flux.mass += dtj * A * totflux[0];
pj->conserved.flux.momentum[0] += dtj * A * totflux[1];
pj->conserved.flux.momentum[1] += dtj * A * totflux[2];
diff --git a/src/hydro/Shadowswift/hydro_part.h b/src/hydro/Shadowswift/hydro_part.h
index 7a5724c3300da8950022dad76f57ae3fa4b9e807..36743bb7ad68123eea4994f1acd3a5a36f283f19 100644
--- a/src/hydro/Shadowswift/hydro_part.h
+++ b/src/hydro/Shadowswift/hydro_part.h
@@ -37,7 +37,7 @@ struct xpart {
/* Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
-} __attribute__((aligned(xpart_align)));
+} SWIFT_STRUCT_ALIGN;
/* Data of a single particle. */
struct part {
@@ -54,12 +54,6 @@ struct part {
/* Particle cutoff radius. */
float h;
- /* Particle time of beginning of time-step. */
- int ti_begin;
-
- /* Particle time of end of time-step. */
- int ti_end;
-
/* The primitive hydrodynamical variables. */
struct {
@@ -161,6 +155,9 @@ struct part {
/* Physical time step of the particle. */
float dt;
+ /* Active flag. */
+ char active;
+
/* Actual velocity of the particle. */
float v_full[3];