Add a new parameter to the gravity calculation to control the frequency of the...

Add a new parameter to the gravity calculation to control the frequency of the tree rebuild using a fraction of the number of gpart updates.

examples/parameter_example.yml

@@ -37,9 +37,10 @@ SPH:
 # Parameters for the self-gravity scheme
 Gravity:
   eta:          0.025               # Constant dimensionless multiplier for time integration.
-  theta:        0.7                 # Opening angle (Multipole acceptance criterion)
+  theta:        0.7                 # Opening angle (Multipole acceptance criterion).
   comoving_softening:     0.0026994 # Comoving softening length (in internal units).
   max_physical_softening: 0.0007    # Physical softening length (in internal units).
+  rebuild_frequency:      0.01      # (Optional) Frequency of the gravity-tree rebuild in units of the number of g-particles (this is the default value).
   a_smooth:     1.25                # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
   r_cut_max:    4.5                 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
   r_cut_min:    0.1                 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).

src/engine.c

@@ -3854,6 +3854,11 @@ void engine_rebuild(struct engine *e, int clean_smoothing_length_values) {
   /* Print the status of the system */
   if (e->verbose) engine_print_task_counts(e);
 
+  /* Clear the counters of updates since the last rebuild */
+  e->updates_since_rebuild = 0;
+  e->g_updates_since_rebuild = 0;
+  e->s_updates_since_rebuild = 0;
+
   /* Flag that a rebuild has taken place */
   e->step_props |= engine_step_prop_rebuild;
 
@@ -4543,6 +4548,11 @@ void engine_step(struct engine *e) {
     fflush(e->file_timesteps);
   }
 
+  /* Update the counters */
+  e->updates_since_rebuild += e->updates;
+  e->g_updates_since_rebuild += e->g_updates;
+  e->s_updates_since_rebuild += e->s_updates;
+
   /* Move forward in time */
   e->ti_old = e->ti_current;
   e->ti_current = e->ti_end_min;
@@ -4623,9 +4633,10 @@ void engine_step(struct engine *e) {
     gravity_exact_force_check(e->s, e, 1e-1);
 #endif
 
-  /* Let's trigger a non-SPH rebuild every-so-often for good measure */
-  if (!(e->policy & engine_policy_hydro) &&  // MATTHIEU improve this
-      (e->policy & engine_policy_self_gravity) && e->step % 20 == 0)
+  /* Trigger a tree-rebuild if we passed the frequency threshold */
+  if ((e->policy & engine_policy_self_gravity) &&
+      (e->g_updates_since_rebuild >
+       e->total_nr_gparts * e->gravity_properties->rebuild_frequency))
     e->forcerebuild = 1;
 
   /* Collect the values of rebuild from all nodes and recover the (integer)

src/engine.h

@@ -188,6 +188,11 @@ struct engine {
   /* Number of particles updated in the previous step */
   size_t updates, g_updates, s_updates;
 
+  /* Number of updates since the last rebuild */
+  size_t updates_since_rebuild;
+  size_t g_updates_since_rebuild;
+  size_t s_updates_since_rebuild;
+
   /* Properties of the previous step */
   int step_props;
 

src/gravity_properties.c

@@ -35,11 +35,20 @@
 #define gravity_props_default_a_smooth 1.25f
 #define gravity_props_default_r_cut_max 4.5f
 #define gravity_props_default_r_cut_min 0.1f
+#define gravity_props_default_rebuild_frequency 0.01f
 
 void gravity_props_init(struct gravity_props *p,
                         const struct swift_params *params,
                         const struct cosmology *cosmo) {
 
+  /* Tree updates */
+  p->rebuild_frequency =
+      parser_get_opt_param_float(params, "Gravity:rebuild_frequency",
+                                 gravity_props_default_rebuild_frequency);
+
+  if (p->rebuild_frequency < 0.f || p->rebuild_frequency > 1.f)
+    error("Invalid tree rebuild frequency. Must be in [0., 1.]");
+
   /* Tree-PM parameters */
   p->a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth",
                                            gravity_props_default_a_smooth);
@@ -116,12 +125,16 @@ void gravity_props_print(const struct gravity_props *p) {
 
   message("Self-gravity tree cut-off: r_cut_max=%f", p->r_cut_max);
   message("Self-gravity truncation cut-off: r_cut_min=%f", p->r_cut_min);
+
+  message("Self-gravity tree update frequency: f=%f", p->rebuild_frequency);
 }
 
 #if defined(HAVE_HDF5)
 void gravity_props_print_snapshot(hid_t h_grpgrav,
                                   const struct gravity_props *p) {
 
+  io_write_attribute_f(h_grpgrav, "Tree update frequency",
+                       p->rebuild_frequency);
   io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta);
   io_write_attribute_f(
       h_grpgrav, "Comoving softening length",

src/gravity_properties.h

@@ -36,6 +36,9 @@
  */
 struct gravity_props {
 
+  /*! Frequency of tree-rebuild in units of #gpart updates. */
+  float rebuild_frequency;
+
   /*! Mesh smoothing scale in units of top-level cell size */
   float a_smooth;