From b9028bc7a3907db62a9f75d63832bc19c3bb6a40 Mon Sep 17 00:00:00 2001
From: Darwin <darwin.roduit@unige.ch>
Date: Sat, 15 Mar 2025 21:59:27 +0000
Subject: [PATCH] GEAR sink: Small fixes

---
 .../SinkParticles/HomogeneousBoxSinkParticles/makeIC.py  | 3 ++-
 src/sink/GEAR/sink.h                                     | 9 ++++++++-
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/makeIC.py b/examples/SinkParticles/HomogeneousBoxSinkParticles/makeIC.py
index 6ad3be507b..80ee002e2d 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 f2b5c00099..47baafcffc 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;
-- 
GitLab