diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h index 459de24ef3c5e9140fd136155ba55d2364795fb8..f438425b540b13a9d637a942b4c12bdcf4cd8c08 100644 --- a/src/chemistry/EAGLE/chemistry.h +++ b/src/chemistry/EAGLE/chemistry.h @@ -78,6 +78,22 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density( struct part* restrict p, const struct chemistry_global_data* cd, const struct cosmology* cosmo) {} +/** + * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. + * + * @param p The particle to act upon + * @param xp The extended particle data to act upon + * @param cd #chemistry_global_data containing chemistry informations. + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void +chemistry_part_has_no_neighbours(struct part* restrict p, + struct xpart* restrict xp, + const struct chemistry_global_data* cd, + const struct cosmology* cosmo) { + error("Needs implementing!"); +} + /** * @brief Sets the chemistry properties of the (x-)particles to a valid start * state. diff --git a/src/chemistry/GEAR/chemistry.h b/src/chemistry/GEAR/chemistry.h index a51051ca3ae45476986d39c868a9fc71bf7f9ae5..b138e66a844b4dca98ba5bd7d46643f0683385b5 100644 --- a/src/chemistry/GEAR/chemistry.h +++ b/src/chemistry/GEAR/chemistry.h @@ -134,6 +134,22 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density( } } +/** + * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. + * + * @param p The particle to act upon + * @param xp The extended particle data to act upon + * @param cd #chemistry_global_data containing chemistry informations. + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void +chemistry_part_has_no_neighbours(struct part* restrict p, + struct xpart* restrict xp, + const struct chemistry_global_data* cd, + const struct cosmology* cosmo) { + error("Needs implementing!"); +} + /** * @brief Sets the chemistry properties of the (x-)particles to a valid start * state. diff --git a/src/chemistry/none/chemistry.h b/src/chemistry/none/chemistry.h index 3ca51660ddfeead2b7ad0010979b719e59c4934e..6c2e7a21a81ce7b96723904f299378a2e5348a81 100644 --- a/src/chemistry/none/chemistry.h +++ b/src/chemistry/none/chemistry.h @@ -86,6 +86,20 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density( struct part* restrict p, const struct chemistry_global_data* cd, const struct cosmology* cosmo) {} +/** + * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. + * + * @param p The particle to act upon + * @param xp The extended particle data to act upon + * @param cd #chemistry_global_data containing chemistry informations. + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void +chemistry_part_has_no_neighbours(struct part* restrict p, + struct xpart* restrict xp, + const struct chemistry_global_data* cd, + const struct cosmology* cosmo) {} + /** * @brief Sets the chemistry properties of the (x-)particles to a valid start * state. diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h index b1a999b63143437cab8518cfdd96885533d7401e..2c3a9c46f0500fb20aa3cfa2e5feb682b3dcec63 100644 --- a/src/hydro/Default/hydro.h +++ b/src/hydro/Default/hydro.h @@ -339,7 +339,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( /* Re-set problematic values */ p->rho = p->mass * kernel_root * h_inv_dim; - p->density.wcount = kernel_root * kernel_norm * h_inv_dim; + p->density.wcount = kernel_root * h_inv_dim; p->rho_dh = 0.f; p->density.wcount_dh = 0.f; p->density.div_v = 0.f; @@ -423,7 +423,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( /* Reset the time derivatives. */ p->force.u_dt = 0.0f; p->force.h_dt = 0.0f; - p->force.v_sig = 0.0f; + p->force.v_sig = p->force.soundspeed; } /** diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index bc06a24e2a8245556a1042f2459273b8d750489e..26e3bf97dd1924abbe7380d1eaadce75213344df 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -349,7 +349,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( /* Re-set problematic values */ p->rho = p->mass * kernel_root * h_inv_dim; - p->density.wcount = kernel_root * kernel_norm * h_inv_dim; + p->density.wcount = kernel_root * h_inv_dim; p->density.rho_dh = 0.f; p->density.wcount_dh = 0.f; p->density.div_v = 0.f; diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h index 9c4be6af359d7236e483b712065b357c6ed35402..ae0fe13bc139175c62e8b8ee9ab8c78358aff4e8 100644 --- a/src/hydro/GizmoMFM/hydro.h +++ b/src/hydro/GizmoMFM/hydro.h @@ -398,7 +398,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ /* Re-set problematic values */ - p->density.wcount = kernel_root * kernel_norm * h_inv_dim; + p->density.wcount = kernel_root * h_inv_dim; p->density.wcount_dh = 0.f; p->geometry.volume = 1.0f; p->geometry.matrix_E[0][0] = 1.0f; diff --git a/src/hydro/GizmoMFV/hydro.h b/src/hydro/GizmoMFV/hydro.h index 1d5abeaaf63d88b02817a691f160c537d0b1915b..e94a87408e561773d3dc95e2443538e5822510d8 100644 --- a/src/hydro/GizmoMFV/hydro.h +++ b/src/hydro/GizmoMFV/hydro.h @@ -398,7 +398,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ /* Re-set problematic values */ - p->density.wcount = kernel_root * kernel_norm * h_inv_dim; + p->density.wcount = kernel_root * h_inv_dim; p->density.wcount_dh = 0.f; p->geometry.volume = 1.0f; p->geometry.matrix_E[0][0] = 1.0f; diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h index 3f9d99683bde4ed6db64d8aaa5b111e2f67f0969..812f8ad72de55ad7990ee6ef88223a401780bc4b 100644 --- a/src/hydro/Minimal/hydro.h +++ b/src/hydro/Minimal/hydro.h @@ -367,7 +367,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( /* Re-set problematic values */ p->rho = p->mass * kernel_root * h_inv_dim; - p->density.wcount = kernel_root * kernel_norm * h_inv_dim; + p->density.wcount = kernel_root * h_inv_dim; p->density.rho_dh = 0.f; p->density.wcount_dh = 0.f; } @@ -426,7 +426,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( /* Reset the time derivatives. */ p->u_dt = 0.0f; p->force.h_dt = 0.0f; - p->force.v_sig = 0.0f; + p->force.v_sig = p->force.soundspeed; } /** diff --git a/src/hydro/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h index 541029ee06dd2799443fc89b688d7baca3fae0f8..73ffc26b8acf687a5445591ddccd72ea8e8fa8ae 100644 --- a/src/hydro/Minimal/hydro_debug.h +++ b/src/hydro/Minimal/hydro_debug.h @@ -36,16 +36,17 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle( const struct part* p, const struct xpart* xp) { printf( - "x=[%.3e,%.3e,%.3e], " - "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], " - "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n" - "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, " - "time_bin=%d\n", + "\n " + "x=[%.6g, %.6g, %.6g], v=[%.3g, %.3g, %.3g], \n " + "v_full=[%.3g, %.3g, %.3g], a=[%.3g, %.3g, %.3g], \n " + "m=%.3g, u=%.3g, du/dt=%.3g, P=%.3g, c_s=%.3g, \n " + "v_sig=%.3g, h=%.3g, dh/dt=%.3g, wcount=%.3g, rho=%.3g, \n " + "dh_drho=%.3g, time_bin=%d \n", p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0], xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2], - p->u, p->u_dt, p->force.v_sig, hydro_get_comoving_pressure(p), p->h, - p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho, - p->time_bin); + p->mass, p->u, p->u_dt, hydro_get_comoving_pressure(p), + p->force.soundspeed, p->force.v_sig, p->h, p->force.h_dt, + p->density.wcount, p->rho, p->density.rho_dh, p->time_bin); } #endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */ diff --git a/src/hydro/MinimalMultiMat/hydro.h b/src/hydro/MinimalMultiMat/hydro.h index 5383ffda8fe67a591691766e4150e75a9dbd4cb0..cfad6b2b2b389da9f423540cb30f1df4cebc5416 100644 --- a/src/hydro/MinimalMultiMat/hydro.h +++ b/src/hydro/MinimalMultiMat/hydro.h @@ -368,7 +368,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( /* Re-set problematic values */ p->rho = p->mass * kernel_root * h_inv_dim; - p->density.wcount = kernel_root * kernel_norm * h_inv_dim; + p->density.wcount = kernel_root * h_inv_dim; p->density.rho_dh = 0.f; p->density.wcount_dh = 0.f; } @@ -429,7 +429,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( /* Reset the time derivatives. */ p->u_dt = 0.0f; p->force.h_dt = 0.0f; - p->force.v_sig = 0.0f; + p->force.v_sig = p->force.soundspeed; } /** diff --git a/src/hydro/MinimalMultiMat/hydro_debug.h b/src/hydro/MinimalMultiMat/hydro_debug.h index d8fe73313fbb175faa9547970c6680346dad0a1b..17b624ad0f660152be4ba685905a3c855e1761f8 100644 --- a/src/hydro/MinimalMultiMat/hydro_debug.h +++ b/src/hydro/MinimalMultiMat/hydro_debug.h @@ -38,16 +38,17 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle( const struct part* p, const struct xpart* xp) { printf( - "x=[%.3e,%.3e,%.3e], " - "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], " - "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n" - "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, " - "time_bin=%d, mat_id=%d\n", + "\n " + "x=[%.6g, %.6g, %.6g], v=[%.3g, %.3g, %.3g], \n " + "v_full=[%.3g, %.3g, %.3g], a=[%.3g, %.3g, %.3g], \n " + "m=%.3g, u=%.3g, du/dt=%.3g, P=%.3g, c_s=%.3g, \n " + "v_sig=%.3g, h=%.3g, dh/dt=%.3g, wcount=%.3g, rho=%.3g, \n " + "dh_drho=%.3g, time_bin=%d, mat_id=%d \n", p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0], xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2], - p->u, p->u_dt, p->force.v_sig, hydro_get_comoving_pressure(p), p->h, - p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho, - p->time_bin, p->mat_id); + p->mass, p->u, p->u_dt, hydro_get_comoving_pressure(p), + p->force.soundspeed, p->force.v_sig, p->h, p->force.h_dt, + p->density.wcount, p->rho, p->density.rho_dh, p->time_bin, p->mat_id); } #endif /* SWIFT_MINIMAL_MULTI_MAT_HYDRO_DEBUG_H */ diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h index 1e0d8208e82f7f48691d1df7603a6b02d1471c12..ea086daeeb1e93d7f1476302564fb4182a6fb611 100644 --- a/src/hydro/PressureEnergy/hydro.h +++ b/src/hydro/PressureEnergy/hydro.h @@ -403,7 +403,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( p->rho = p->mass * kernel_root * h_inv_dim; p->pressure_bar = p->mass * p->u * hydro_gamma_minus_one * kernel_root * h_inv_dim; - p->density.wcount = kernel_root * kernel_norm * h_inv_dim; + p->density.wcount = kernel_root * h_inv_dim; p->density.rho_dh = 0.f; p->density.wcount_dh = 0.f; p->density.pressure_bar_dh = 0.f; @@ -480,7 +480,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( /* Reset the time derivatives. */ p->u_dt = 0.0f; p->force.h_dt = 0.0f; - p->force.v_sig = 0.0f; + p->force.v_sig = p->force.soundspeed; } /** diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h index 87d46c6d43f0d4f6de6d18f5400b38f0fc4d0f55..e4b7cf06e083638a94526cc1f9e7212cf19dfad4 100644 --- a/src/hydro/PressureEntropy/hydro.h +++ b/src/hydro/PressureEntropy/hydro.h @@ -357,7 +357,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( /* Re-set problematic values */ p->rho = p->mass * kernel_root * h_inv_dim; p->rho_bar = p->mass * kernel_root * h_inv_dim; - p->density.wcount = kernel_root * kernel_norm * h_inv_dim; + p->density.wcount = kernel_root * h_inv_dim; p->density.rho_dh = 0.f; p->density.wcount_dh = 0.f; p->density.pressure_dh = 0.f; @@ -441,7 +441,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( p->force.h_dt = 0.0f; /* Reset maximal signal velocity */ - p->force.v_sig = 0.0f; + p->force.v_sig = p->force.soundspeed; } /** diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h index 36078798cdd1ac68a456134fd7887408752f18c9..025779f17496e7bf30fdf12353c4381c7d6292ce 100644 --- a/src/hydro/Shadowswift/hydro.h +++ b/src/hydro/Shadowswift/hydro.h @@ -257,7 +257,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ /* Re-set problematic values */ - p->density.wcount = kernel_root * kernel_norm * h_inv_dim; + p->density.wcount = kernel_root * h_inv_dim; p->density.wcount_dh = 0.f; } diff --git a/src/runner.c b/src/runner.c index 4409411cf7fac01d2358f182f14a2d00dc02cdc6..4ab2075f34aa1a4bdaea70e29c2d4894b948bfb3 100644 --- a/src/runner.c +++ b/src/runner.c @@ -710,9 +710,13 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { const float h_old_dim = pow_dimension(h_old); const float h_old_dim_minus_one = pow_dimension_minus_one(h_old); float h_new; + int has_no_neighbours = 0; if (p->density.wcount == 0.f) { /* No neighbours case */ + /* Flag that there were no neighbours */ + has_no_neighbours = 1; + /* Double h and try again */ h_new = 2.f * h_old; } else { @@ -729,7 +733,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { p->density.wcount_dh * h_old_dim + hydro_dimension * p->density.wcount * h_old_dim_minus_one; - h_new = h_old - f / f_prime; + /* Avoid floating point exception from f_prime = 0 */ + h_new = h_old - f / (f_prime + FLT_MIN); #ifdef SWIFT_DEBUG_CHECKS if ((f > 0.f && h_new > h_old) || (f < 0.f && h_new < h_old)) @@ -768,8 +773,10 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { p->h = hydro_h_max; /* Do some damage control if no neighbours at all were found */ - if (p->density.wcount == kernel_root * kernel_norm) + if (has_no_neighbours) { hydro_part_has_no_neighbours(p, xp, cosmo); + chemistry_part_has_no_neighbours(p, xp, chemistry, cosmo); + } } }