diff --git a/src/const.h b/src/const.h index 92ddce19dcfb4e4c8af583c28a605c821d5bdee4..1239d8cd6894aa678f706bc2f67df764e259d533 100644 --- a/src/const.h +++ b/src/const.h @@ -57,7 +57,7 @@ /* Options to control the movement of particles for GIZMO_SPH. */ /* This option disables particle movement */ -#define GIZMO_FIX_PARTICLES +//#define GIZMO_FIX_PARTICLES /* Source terms */ #define SOURCETERMS_NONE diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h index cb2978e8f5447c109db3e6cdbc9e022eb3ef28ba..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) @@ -79,10 +80,20 @@ __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) { + 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.; @@ -189,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; + } } /** @@ -273,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 @@ -363,13 +354,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( /* Add normalization to h_dt. */ p->force.h_dt *= p->h * hydro_dimension_inv; -#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; -#else /* 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 @@ -396,6 +380,17 @@ __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 611f0449a5951f9d5ed2667b65be043dcedcf8f3..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 {