diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/makeIC.py b/examples/SinkParticles/HomogeneousBoxSinkParticles/makeIC.py index 6ad3be507b4fd12d8bef3c0af4a1b674cb950a27..80ee002e2de865068e0a1190a9814d7c12f58d48 100755 --- a/examples/SinkParticles/HomogeneousBoxSinkParticles/makeIC.py +++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/makeIC.py @@ -201,7 +201,6 @@ N_sink = opt.n_sink m_sink = opt.sink_mass * units.M_sun m_sink = m_sink.to(UnitMass).value # Convert the sink mass to internal units - if N_sink == 1: pos_sink = np.reshape(opt.sinks_vel, (N_sink, 3)) else: @@ -215,6 +214,7 @@ else: np.zeros([N_sink, 3]) mass_sink = np.ones(N_sink) * m_sink +h_sink = np.ones(N_sink) * 3 * L / (N + N_sink) ** (1 / 3.0) ids_sink = np.arange(N, N + N_sink) ##################### @@ -263,5 +263,6 @@ grp = fileOutput.create_group("/PartType3") grp.create_dataset("Coordinates", data=pos_sink, dtype="d") grp.create_dataset("Velocities", data=vel_sink, dtype="f") grp.create_dataset("Masses", data=mass_sink, dtype="f") +grp.create_dataset("SmoothingLength", data=h_sink, dtype="f") grp.create_dataset("ParticleIDs", data=ids_sink, dtype="L") fileOutput.close() diff --git a/src/sink/GEAR/sink.h b/src/sink/GEAR/sink.h index f2b5c000995e74b7d4a7f404e0170a710de482ea..47baafcffc08b4645ef8679427ea73433743a21a 100644 --- a/src/sink/GEAR/sink.h +++ b/src/sink/GEAR/sink.h @@ -77,7 +77,8 @@ __attribute__((always_inline)) INLINE static float sink_compute_timestep( sink->to_collect.sound_speed_gas * cosmo->a_factor_sound_speed; const double gas_c_phys2 = gas_c_phys * gas_c_phys; const float denominator = sqrtf(gas_c_phys2 + gas_v_norm2); - const float h_min = cosmo->a * min(sink->h, sink->to_collect.minimal_h_gas); + const float h_min = + cosmo->a * kernel_gamma * min(sink->h, sink->to_collect.minimal_h_gas); float dt_cfl = 0.0; /* This case can happen if the sink is just born. */ @@ -165,6 +166,12 @@ __attribute__((always_inline)) INLINE static void sink_first_init_sink( struct sink* sp, const struct sink_props* sink_props, const struct engine* e) { + /* Set the smoothing length if it is fixed. Note that, otherwise, the + smoothing lengths were read from the ICs. */ + if (sink_props->use_fixed_r_cut) { + sp->h = sink_props->cut_off_radius / kernel_gamma; + } + sp->time_bin = 0; sp->number_of_gas_swallows = 0;