Future-proof version of dokick that also gets time-step info from self-gravity.

examples/ExternalPointMass/externalPointMass.yml

@@ -28,7 +28,7 @@ SPH:
   delta_neighbours:      1.       # The tolerance for the targetted number of neighbours.
   CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
   max_ghost_iterations:  30       # Maximal number of iterations allowed to converge towards the smoothing length.
-  max_smoothing_length:  3.       # Maximal smoothing length allowed (in internal units).
+  max_smoothing_length:  10.      # Maximal smoothing length allowed (in internal units).
 
 # Parameters related to the initial conditions
 InitialConditions:

src/gravity/Default/gravity.h

@@ -22,7 +22,8 @@
 #include "potentials.h"
 
 /**
- * @brief Computes the gravity time-step of a given particle.
+ * @brief Computes the gravity time-step of a given particle due to an external
+ *potential.
  *
  * This function only branches towards the potential chosen by the user.
  *
@@ -30,20 +31,40 @@
  * @param phys_const The physical constants in internal units.
  * @param gp Pointer to the g-particle data.
  */
-__attribute__((always_inline)) INLINE static float gravity_compute_timestep(
-    const struct external_potential* potential,
-    const struct phys_const* const phys_const, const struct gpart* const gp) {
+__attribute__((always_inline))
+    INLINE static float gravity_compute_timestep_external(
+        const struct external_potential* potential,
+        const struct phys_const* const phys_const,
+        const struct gpart* const gp) {
 
   float dt = FLT_MAX;
 
 #ifdef EXTERNAL_POTENTIAL_POINTMASS
-  dt =   
-      fminf(dt,  external_gravity_pointmass_timestep(potential, phys_const, gp));
+  dt =
+      fminf(dt, external_gravity_pointmass_timestep(potential, phys_const, gp));
 #endif
 
   return dt;
 }
 
+/**
+ * @brief Computes the gravity time-step of a given particle due to self-gravity
+ *
+ * This function only branches towards the potential chosen by the user.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param gp Pointer to the g-particle data.
+ */
+__attribute__((always_inline))
+    INLINE static float gravity_compute_timestep_self(
+        const struct phys_const* const phys_const,
+        const struct gpart* const gp) {
+
+  float dt = FLT_MAX;
+
+  return dt;
+}
+
 /**
  * @brief Initialises the g-particles for the first time
  *

src/runner.c

@@ -111,6 +111,8 @@ void runner_dograv_external(struct runner *r, struct cell *c, int timer) {
   struct gpart *g, *gparts = c->gparts;
   int i, k, gcount = c->gcount;
   const int ti_current = r->e->ti_current;
+  const struct external_potential *potential = r->e->potential;
+  const struct phys_const *constants = r->e->physical_constants;
 
   TIMER_TIC;
 
@@ -134,7 +136,7 @@ void runner_dograv_external(struct runner *r, struct cell *c, int timer) {
     /* Is this part within the time step? */
     if (g->ti_end <= ti_current) {
 
-      external_gravity(r->e->potential, r->e->physical_constants, g);
+      external_gravity(potential, constants, g);
 
       /* /\* check for energy and angular momentum conservation - begin by */
       /*  * synchronizing velocity*\/ */
@@ -900,6 +902,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
   struct part *const parts = c->parts;
   struct xpart *const xparts = c->xparts;
   struct gpart *const gparts = c->gparts;
+  const struct external_potential *potential = r->e->potential;
+  const struct phys_const *constants = r->e->physical_constants;
   const int is_fixdt =
       (r->e->policy & engine_policy_fixdt) == engine_policy_fixdt;
 
@@ -941,8 +945,12 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
         } else {
 
           /* Compute the next timestep (gravity condition) */
-          float new_dt = gravity_compute_timestep(r->e->potential,
-                                                  r->e->physical_constants, gp);
+          const float new_dt_external =
+              gravity_compute_timestep_external(potential, constants, gp);
+          const float new_dt_self =
+              gravity_compute_timestep_self(constants, gp);
+
+          float new_dt = fminf(new_dt_external, new_dt_self);
 
           /* Limit timestep within the allowed range */
           new_dt = fminf(new_dt, global_dt_max);
@@ -1026,10 +1034,17 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
 
           /* Compute the next timestep (gravity condition) */
           float new_dt_grav = FLT_MAX;
-          if (p->gpart != NULL)
-            new_dt_grav = gravity_compute_timestep(
-                r->e->potential, r->e->physical_constants, p->gpart);
+          if (p->gpart != NULL) {
+
+            const float new_dt_external = gravity_compute_timestep_external(
+                potential, constants, p->gpart);
+            const float new_dt_self =
+                gravity_compute_timestep_self(constants, p->gpart);
+
+            new_dt_grav = fminf(new_dt_external, new_dt_self);
+          }
 
+          /* Final time-step is minimum of hydro and gravity */
           float new_dt = fminf(new_dt_hydro, new_dt_grav);
 
           /* Limit change in h */

src/space.c

@@ -177,10 +177,6 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
       }
       s->h_max = h_max;
     }
-  } else {
-    /* It would be nice to replace this with something more physical or
-     * meaningful */
-    h_max = s->dim[0] / 16.f;
   }
 
 /* If we are running in parallel, make sure everybody agrees on