diff --git a/examples/EAGLE_100/eagle_100.yml b/examples/EAGLE_100/eagle_100.yml index 3bcda091bdac5b740f3568de9c0814cc84c3b846..5275f709c710d2da1d5396e98f1d4d918e482c6d 100644 --- a/examples/EAGLE_100/eagle_100.yml +++ b/examples/EAGLE_100/eagle_100.yml @@ -49,6 +49,7 @@ Gravity: # Parameters for the hydrodynamics scheme SPH: resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). + h_min_ratio: 0.1 # Minimal smoothing in units of softening. CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. minimal_temperature: 100 # (internal units) diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml index 367c5f769775d9293c0ec7b0bf670a2b22dcb328..90070f6a342d60d48d8fd2c794014f9d638030be 100644 --- a/examples/EAGLE_12/eagle_12.yml +++ b/examples/EAGLE_12/eagle_12.yml @@ -47,6 +47,7 @@ Gravity: # Parameters for the hydrodynamics scheme SPH: resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). + h_min_ratio: 0.1 # Minimal smoothing in units of softening. CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. minimal_temperature: 100 # (internal units) diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml index cab0dbcd5efc0528ddc65a6dde1e5c2d7cb6b9a9..3d7d5de4fe2c98a04c986840bc2804edaf780c17 100644 --- a/examples/EAGLE_25/eagle_25.yml +++ b/examples/EAGLE_25/eagle_25.yml @@ -54,6 +54,7 @@ Gravity: # Parameters for the hydrodynamics scheme SPH: resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). + h_min_ratio: 0.1 # Minimal smoothing in units of softening. CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. minimal_temperature: 100 # (internal units) diff --git a/examples/EAGLE_50/eagle_50.yml b/examples/EAGLE_50/eagle_50.yml index 35f7c85993d6d52189eaa1da35821fc5fa103c56..8d622998c5ea1c754e877643e42ee3d97bfa720e 100644 --- a/examples/EAGLE_50/eagle_50.yml +++ b/examples/EAGLE_50/eagle_50.yml @@ -49,6 +49,7 @@ Gravity: # Parameters for the hydrodynamics scheme SPH: resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). + h_min_ratio: 0.1 # Minimal smoothing in units of softening. CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. minimal_temperature: 100 # (internal units) diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml index e80fac8167a832c17cd10e1d2ae7cd854f314d17..fb885fa44f119bb3e296947997f71581d5ca0f48 100644 --- a/examples/EAGLE_6/eagle_6.yml +++ b/examples/EAGLE_6/eagle_6.yml @@ -58,6 +58,7 @@ Gravity: # Parameters for the hydrodynamics scheme SPH: resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). + h_min_ratio: 0.1 # Minimal smoothing in units of softening. CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. minimal_temperature: 100 # (internal units) diff --git a/examples/IsolatedGalaxy_starformation/isolated_galaxy.yml b/examples/IsolatedGalaxy_starformation/isolated_galaxy.yml index d518732056f508cbf896dcd5f258663522a91eea..2a885f5770f88c3ee6dd75e194d0fa713160e39e 100644 --- a/examples/IsolatedGalaxy_starformation/isolated_galaxy.yml +++ b/examples/IsolatedGalaxy_starformation/isolated_galaxy.yml @@ -11,13 +11,13 @@ Gravity: mesh_side_length: 32 # Number of cells along each axis for the periodic gravity mesh. eta: 0.025 # Constant dimensionless multiplier for time integration. theta: 0.7 # Opening angle (Multipole acceptance criterion). - comoving_softening: 0.01 # Comoving softening length (in internal units). - max_physical_softening: 0.01 # Physical softening length (in internal units). + comoving_softening: 0.01 # Comoving softening length (in internal units). + max_physical_softening: 0.01 # Physical softening length (in internal units). # Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.) TimeIntegration: time_begin: 0. # The starting time of the simulation (in internal units). - time_end: 0.1 # 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-9 # 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). @@ -35,20 +35,17 @@ Statistics: # Parameters related to the initial conditions InitialConditions: - file_name: fid.hdf5 # The file to read + file_name: fid.hdf5 # The file to read periodic: 0 # Are we running with periodic ICs? # Parameters for the hydrodynamics scheme SPH: resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. - h_tolerance: 1e-4 # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths. - h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified. - max_volume_change: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step. - max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length. - minimal_temperature: 100 # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0. - H_ionization_temperature: 1e4 # (Optional) Temperature of the transition from neutral to ionized Hydrogen for primoridal gas. + h_min_ratio: 0.1 # Minimal smoothing in units of softening. + h_max: 10. +# Standard EAGLE cooling options EAGLECooling: dir_name: ./coolingtables/ # Location of the Wiersma+08 cooling tables H_reion_z: 11.5 # Redshift of Hydrogen re-ionization diff --git a/examples/SantaBarbara/santa_barbara.yml b/examples/SantaBarbara/santa_barbara.yml index eb20030e00c36646b4283012eab5ecae0e286680..aa478a29023af1945bdd8cf34ace1c421f249fc1 100644 --- a/examples/SantaBarbara/santa_barbara.yml +++ b/examples/SantaBarbara/santa_barbara.yml @@ -46,6 +46,7 @@ Gravity: # Parameters of the hydro scheme SPH: resolution_eta: 1.2348 # "48 Ngb" with the cubic spline kernel + h_min_ratio: 0.1 CFL_condition: 0.1 initial_temperature: 1200. # (1 + z_ini)^2 * 2.72K minimal_temperature: 100. diff --git a/examples/SantaBarbara_low/santa_barbara.yml b/examples/SantaBarbara_low/santa_barbara.yml index 0e3c66b1a1c1e04bf6fad30e806327b83f03737e..35a6e2762f91707ed43bd5f0107c3dd8a53e12e4 100644 --- a/examples/SantaBarbara_low/santa_barbara.yml +++ b/examples/SantaBarbara_low/santa_barbara.yml @@ -47,6 +47,7 @@ Gravity: SPH: resolution_eta: 1.2348 # "48 Ngb" with the cubic spline kernel CFL_condition: 0.1 + h_min_ratio: 0.1 initial_temperature: 1200. # (1 + z_ini)^2 * 2.72K minimal_temperature: 100. diff --git a/examples/SmallCosmoVolume/small_cosmo_volume.yml b/examples/SmallCosmoVolume/small_cosmo_volume.yml index 32846227c93e7d747fade2dc306b4c5c5246c14f..e0d633079e941ade161b7e2fde0fbc063cbac254 100644 --- a/examples/SmallCosmoVolume/small_cosmo_volume.yml +++ b/examples/SmallCosmoVolume/small_cosmo_volume.yml @@ -30,6 +30,7 @@ Gravity: # Parameters of the hydro scheme SPH: resolution_eta: 1.2348 # "48 Ngb" with the cubic spline kernel + h_min_ratio: 0.1 CFL_condition: 0.1 initial_temperature: 7075. # (1 + z_ini)^2 * 2.72K minimal_temperature: 100. diff --git a/examples/SmallCosmoVolume_VELOCIraptor/small_cosmo_volume.yml b/examples/SmallCosmoVolume_VELOCIraptor/small_cosmo_volume.yml index d6b9a78fe3c2a891492affbdea9787d62916d3ed..15007ca8d39a328166a208a2d4cbbc6aea580009 100644 --- a/examples/SmallCosmoVolume_VELOCIraptor/small_cosmo_volume.yml +++ b/examples/SmallCosmoVolume_VELOCIraptor/small_cosmo_volume.yml @@ -30,6 +30,7 @@ Gravity: # Parameters of the hydro scheme SPH: resolution_eta: 1.2348 # "48 Ngb" with the cubic spline kernel + h_min_ratio: 0.1 CFL_condition: 0.1 initial_temperature: 7075. # (1 + z_ini)^2 * 2.72K minimal_temperature: 100. diff --git a/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml b/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml index bb4546659303a05aef1696fe4f2195031d9c035c..60046cadc3bd3af324b5f967c76315cf5a27fe52 100644 --- a/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml +++ b/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml @@ -30,6 +30,7 @@ Gravity: # Parameters of the hydro scheme SPH: resolution_eta: 1.2348 # "48 Ngb" with the cubic spline kernel + h_min_ratio: 0.1 CFL_condition: 0.1 initial_temperature: 7075. # (1 + z_ini)^2 * 2.72K minimal_temperature: 100. diff --git a/examples/ZeldovichPancake_3D/zeldovichPancake.yml b/examples/ZeldovichPancake_3D/zeldovichPancake.yml index a1d2342b56d6816c3cfbe7da70220ab244104fbd..6a7c5166635b7fa0ed5f69c41461d867c3b254ad 100644 --- a/examples/ZeldovichPancake_3D/zeldovichPancake.yml +++ b/examples/ZeldovichPancake_3D/zeldovichPancake.yml @@ -35,6 +35,7 @@ Statistics: SPH: resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. + h_min_ratio: 0.1 # Parameters related to the initial conditions InitialConditions: diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index dae922378b24614211e38cef22548cf0404ee196..39abf52c571bd65fb9904d7d66ffa662df8e0eec 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -28,6 +28,7 @@ SPH: CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. h_tolerance: 1e-4 # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths. h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified. + h_min_ratio: 0. # (Optional) Minimal allowed smoothing length in units of the softening. Defaults to 0 if unspecified. max_volume_change: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step. max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length. initial_temperature: 0 # (Optional) Initial temperature (in internal units) to set the gas particles at start-up. Value is ignored if set to 0. diff --git a/src/cell.c b/src/cell.c index bd1022f1fa23b5911c4056b602008601fa36ce68..4f44abd924be97e54e52fffda8559f7a0dcfe68a 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3642,6 +3642,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) { const double dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]}; const int with_cosmology = (e->policy & engine_policy_cosmology); const float hydro_h_max = e->hydro_properties->h_max; + const float hydro_h_min = e->hydro_properties->h_min; const integertime_t ti_old_part = c->hydro.ti_old_part; const integertime_t ti_current = e->ti_current; struct part *const parts = c->hydro.parts; @@ -3776,6 +3777,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) { /* Limit h to within the allowed range */ p->h = min(p->h, hydro_h_max); + p->h = max(p->h, hydro_h_min); /* Compute (square of) motion since last cell construction */ const float dx2 = xp->x_diff[0] * xp->x_diff[0] + diff --git a/src/engine.c b/src/engine.c index f32bc0dbe7e2346f9d15984575f1128a0b589d21..4a34112ca4b758768d3171fc49546a7f4d1e49de 100644 --- a/src/engine.c +++ b/src/engine.c @@ -2792,7 +2792,12 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs, /* Update the softening lengths */ if (e->policy & engine_policy_self_gravity) - gravity_update(e->gravity_properties, e->cosmology); + gravity_props_update(e->gravity_properties, e->cosmology); + + /* Udpate the hydro properties */ + if (e->policy & engine_policy_hydro) + hydro_props_update(e->hydro_properties, e->gravity_properties, + e->cosmology); /* Start by setting the particles in a good state */ if (e->nodeID == 0) message("Setting particles to a valid state..."); @@ -3085,7 +3090,12 @@ void engine_step(struct engine *e) { /* Update the softening lengths */ if (e->policy & engine_policy_self_gravity) - gravity_update(e->gravity_properties, e->cosmology); + gravity_props_update(e->gravity_properties, e->cosmology); + + /* Udpate the hydro properties */ + if (e->policy & engine_policy_hydro) + hydro_props_update(e->hydro_properties, e->gravity_properties, + e->cosmology); /* Trigger a tree-rebuild if we passed the frequency threshold */ if ((e->policy & engine_policy_self_gravity) && @@ -4088,7 +4098,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params, int policy, int verbose, struct repartition *reparttype, const struct unit_system *internal_units, const struct phys_const *physical_constants, - struct cosmology *cosmo, const struct hydro_props *hydro, + struct cosmology *cosmo, struct hydro_props *hydro, const struct entropy_floor_properties *entropy_floor, struct gravity_props *gravity, const struct stars_props *stars, struct pm_mesh *mesh, diff --git a/src/engine.h b/src/engine.h index 5a9959a6a024b7e0d2ae85f63ce87302ebe9e0a7..938213d8c5b8889266e74e38f6fb2c00baa1d256 100644 --- a/src/engine.h +++ b/src/engine.h @@ -343,7 +343,7 @@ struct engine { struct cosmology *cosmology; /* Properties of the hydro scheme */ - const struct hydro_props *hydro_properties; + struct hydro_props *hydro_properties; /* Properties of the entropy floor */ const struct entropy_floor_properties *entropy_floor; @@ -419,7 +419,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params, int policy, int verbose, struct repartition *reparttype, const struct unit_system *internal_units, const struct phys_const *physical_constants, - struct cosmology *cosmo, const struct hydro_props *hydro, + struct cosmology *cosmo, struct hydro_props *hydro, const struct entropy_floor_properties *entropy_floor, struct gravity_props *gravity, const struct stars_props *stars, struct pm_mesh *mesh, diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h index d446844e8ffc862fd3be0688302ebb3a2efab8fa..6d1270c952100f3a25202fcdb22be09f9acaa8d9 100644 --- a/src/gravity/Default/gravity.h +++ b/src/gravity/Default/gravity.h @@ -40,7 +40,7 @@ __attribute__((always_inline)) INLINE static float gravity_get_mass( } /** - * @brief Returns the softening of a particle + * @brief Returns the current co-moving softening of a particle * * @param gp The particle of interest * @param grav_props The global gravity properties. diff --git a/src/gravity/Potential/gravity.h b/src/gravity/Potential/gravity.h index 10628dcea91786b0c7483134b8f7f844d6359e49..7d38e9126f1a313b169092b080ebc312c4bbe1bc 100644 --- a/src/gravity/Potential/gravity.h +++ b/src/gravity/Potential/gravity.h @@ -39,7 +39,7 @@ __attribute__((always_inline)) INLINE static float gravity_get_mass( } /** - * @brief Returns the softening of a particle + * @brief Returns the current co-moving softening of a particle * * @param gp The particle of interest * @param grav_props The global gravity properties. diff --git a/src/gravity_properties.c b/src/gravity_properties.c index e548e3010f3b46065a2510723b5bde97121b4c02..60c0bd05ba8cb2234284154b4383bec07f30756d 100644 --- a/src/gravity_properties.c +++ b/src/gravity_properties.c @@ -99,10 +99,11 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params, } /* Set the softening to the current time */ - gravity_update(p, cosmo); + gravity_props_update(p, cosmo); } -void gravity_update(struct gravity_props *p, const struct cosmology *cosmo) { +void gravity_props_update(struct gravity_props *p, + const struct cosmology *cosmo) { /* Current softening lengths */ double softening; diff --git a/src/gravity_properties.h b/src/gravity_properties.h index 0cabd9958efa2bb23524d03632f90fdd1f1c8306..09c8ef8ffa1d6cc4effa4895614106217a2861a9 100644 --- a/src/gravity_properties.h +++ b/src/gravity_properties.h @@ -73,7 +73,7 @@ struct gravity_props { /*! Maxium physical softening */ double epsilon_max_physical; - /*! Current sftening length */ + /*! Current softening length */ float epsilon_cur; /*! Square of current softening length */ @@ -90,7 +90,8 @@ void gravity_props_print(const struct gravity_props *p); void gravity_props_init(struct gravity_props *p, struct swift_params *params, const struct cosmology *cosmo, int with_cosmology, int periodic); -void gravity_update(struct gravity_props *p, const struct cosmology *cosmo); +void gravity_props_update(struct gravity_props *p, + const struct cosmology *cosmo); #if defined(HAVE_HDF5) void gravity_props_print_snapshot(hid_t h_grpsph, diff --git a/src/hydro_properties.c b/src/hydro_properties.c index e83b90cdd932e2a61baaace72922025cf6e199f6..f14c88bfb5128c1da17590f50698e5d038734b71 100644 --- a/src/hydro_properties.c +++ b/src/hydro_properties.c @@ -30,12 +30,16 @@ #include "dimension.h" #include "equation_of_state.h" #include "error.h" +#include "gravity_properties.h" #include "hydro.h" #include "kernel_hydro.h" +#include "parser.h" +#include "units.h" #define hydro_props_default_max_iterations 30 #define hydro_props_default_volume_change 1.4f #define hydro_props_default_h_max FLT_MAX +#define hydro_props_default_h_min_ratio 0.f #define hydro_props_default_h_tolerance 1e-4 #define hydro_props_default_init_temp 0.f #define hydro_props_default_min_temp 0.f @@ -104,6 +108,13 @@ void hydro_props_init(struct hydro_props *p, p->h_max = parser_get_opt_param_float(params, "SPH:h_max", hydro_props_default_h_max); + /* Minimal smoothing length ratio to softening */ + p->h_min_ratio = parser_get_opt_param_float(params, "SPH:h_min_ratio", + hydro_props_default_h_min_ratio); + + /* Temporarily set the minimal softening to 0. */ + p->h_min = 0.f; + /* Number of iterations to converge h */ p->max_smoothing_iterations = parser_get_opt_param_int( params, "SPH:max_ghost_iterations", hydro_props_default_max_iterations); @@ -317,6 +328,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) { * @param p the struct */ void hydro_props_init_no_hydro(struct hydro_props *p) { + p->eta_neighbours = 1.2348; p->h_tolerance = hydro_props_default_h_tolerance; p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm; @@ -325,6 +337,8 @@ void hydro_props_init_no_hydro(struct hydro_props *p) { (pow_dimension(delta_eta) - pow_dimension(p->eta_neighbours)) * kernel_norm; p->h_max = hydro_props_default_h_max; + p->h_min = 0.f; + p->h_min_ratio = hydro_props_default_h_min_ratio; p->max_smoothing_iterations = hydro_props_default_max_iterations; p->CFL_condition = 0.1; p->log_max_h_change = logf(powf(1.4, hydro_dimension_inv)); @@ -352,6 +366,20 @@ void hydro_props_init_no_hydro(struct hydro_props *p) { p->diffusion.alpha_min = hydro_props_default_diffusion_alpha_min; } +/** + * @brief Update the global properties of the hydro scheme for that time-step. + * + * @param p The properties to update. + * @param gp The properties of the gravity scheme. + * @param cosmo The cosmological model. + */ +void hydro_props_update(struct hydro_props *p, const struct gravity_props *gp, + const struct cosmology *cosmo) { + + /* Update the minimal allowed smoothing length */ + p->h_min = p->h_min_ratio * gp->epsilon_cur; +} + /** * @brief Write a hydro_props struct to the given FILE as a stream of bytes. * diff --git a/src/hydro_properties.h b/src/hydro_properties.h index c63ecd9af777b6420858ee2a732d3012c123b209..afc8a4b87aeb39f6bff1e61ae39edf391c856b1b 100644 --- a/src/hydro_properties.h +++ b/src/hydro_properties.h @@ -32,10 +32,14 @@ #endif /* Local includes. */ -#include "parser.h" -#include "physical_constants.h" #include "restart.h" -#include "units.h" + +/* Forward declarations */ +struct cosmology; +struct swift_params; +struct gravity_props; +struct phys_const; +struct unit_system; /** * @brief Contains all the constants and parameters of the hydro scheme @@ -57,6 +61,12 @@ struct hydro_props { /*! Maximal smoothing length */ float h_max; + /*! Minimal smoothing length expressed as ratio to softening length */ + float h_min_ratio; + + /*! Minimal smoothing length */ + float h_min; + /*! Maximal number of iterations to converge h */ int max_smoothing_iterations; @@ -131,6 +141,9 @@ void hydro_props_init(struct hydro_props *p, const struct unit_system *us, struct swift_params *params); +void hydro_props_update(struct hydro_props *p, const struct gravity_props *gp, + const struct cosmology *cosmo); + #if defined(HAVE_HDF5) void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p); #endif diff --git a/src/runner.c b/src/runner.c index df8a52e2cae3f06fb4e59b1cfbce3cf20160f947..08f0432d81ecaef3fbbd4ffcca494868dad960dd 100644 --- a/src/runner.c +++ b/src/runner.c @@ -1230,6 +1230,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { const struct cosmology *cosmo = e->cosmology; const struct chemistry_global_data *chemistry = e->chemistry; const float hydro_h_max = e->hydro_properties->h_max; + const float hydro_h_min = e->hydro_properties->h_min; const float eps = e->hydro_properties->h_tolerance; const float hydro_eta_dim = pow_dimension(e->hydro_properties->eta_neighbours); @@ -1294,6 +1295,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { const float h_old = p->h; 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; @@ -1330,11 +1332,13 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { #endif /* Skip if h is already h_max and we don't have enough neighbours */ - if ((p->h >= hydro_h_max) && (f < 0.f)) { + if (((p->h >= hydro_h_max) && (f < 0.f)) || + ((p->h <= hydro_h_min) && (f > 0.f))) { /* We have a particle whose smoothing length is already set (wants - * to be larger but has already hit the maximum). So, just tidy up - * as if the smoothing length had converged correctly */ + * to be larger but has already hit the maximum OR wants to be smaller + * but has already reached the minimum). So, just tidy up as if the + * smoothing length had converged correctly */ #ifdef EXTRA_HYDRO_LOOP @@ -1436,8 +1440,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { p->h = h_new; } - /* If below the absolute maximum, try again */ - if (p->h < hydro_h_max) { + /* If within the allowed range, try again */ + if (p->h < hydro_h_max && p->h > hydro_h_min) { /* Flag for another round of fun */ pid[redo] = pid[i]; @@ -1453,7 +1457,18 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { /* Off we go ! */ continue; - } else { + } else if (p->h <= hydro_h_min) { + + /* Ok, this particle is a lost cause... */ + p->h = hydro_h_min; + + /* Do some damage control if no neighbours at all were found */ + if (has_no_neighbours) { + hydro_part_has_no_neighbours(p, xp, cosmo); + chemistry_part_has_no_neighbours(p, xp, chemistry, cosmo); + } + + } else if (p->h >= hydro_h_max) { /* Ok, this particle is a lost cause... */ p->h = hydro_h_max; diff --git a/src/space.c b/src/space.c index d930bcdcabada454c86719b66eb029f4e231d6b3..897554225e745819f9f9b9b1fc9c1e35d988713a 100644 --- a/src/space.c +++ b/src/space.c @@ -3351,6 +3351,10 @@ void space_first_init_parts_mapper(void *restrict map_data, int count, const struct hydro_props *hydro_props = s->e->hydro_properties; const float u_init = hydro_props->initial_internal_energy; + const float hydro_h_min_ratio = e->hydro_properties->h_min_ratio; + + const struct gravity_props *grav_props = s->e->gravity_properties; + const int with_gravity = e->policy & engine_policy_self_gravity; const struct chemistry_global_data *chemistry = e->chemistry; const struct cooling_function_data *cool_func = e->cooling_func; @@ -3360,6 +3364,12 @@ void space_first_init_parts_mapper(void *restrict map_data, int count, if (p[k].h <= 0.) error("Invalid value of smoothing length for part %lld h=%e", p[k].id, p[k].h); + + if (with_gravity) { + const struct gpart *gp = p[k].gpart; + const float softening = gravity_get_softening(gp, grav_props); + p->h = max(p->h, softening * hydro_h_min_ratio); + } } /* Convert velocities to internal units */