From dad4582b5f092df16c112251ef3a3fd717871b20 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Fri, 24 Apr 2020 14:41:25 +0200
Subject: [PATCH] Use the geometric MAC on the pre-0 step. Then switch to the
 new MAC if the user asked for it.

---
 examples/EAGLE_low_z/EAGLE_6/eagle_6.yml |  4 ++-
 examples/parameter_example.yml           | 32 ++++++++++++------------
 src/engine.c                             | 25 ++++++++++--------
 src/gravity_properties.c                 | 15 ++++++-----
 src/gravity_properties.h                 |  2 +-
 5 files changed, 43 insertions(+), 35 deletions(-)

diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index 2e7de929fd..baa5f27682 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -54,7 +54,9 @@ Statistics:
 # Parameters for the self-gravity scheme
 Gravity:
   eta:                    0.025    # Constant dimensionless multiplier for time integration.
-  theta:                  0.7      # Opening angle (Multipole acceptance criterion)
+  MAC:                    adaptive
+  epsilon_fmm:            0.01
+  theta_cr:               0.7      # Opening angle (Multipole acceptance criterion)
   mesh_side_length:       16
   dithering:              0
   comoving_DM_softening:         0.0026994 # Comoving DM softening length (in internal units).
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 8310153bf6..af4edf797b 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -62,22 +62,22 @@ Stars:
 
 # Parameters for the self-gravity scheme
 Gravity:
-  mesh_side_length:               128       # Number of cells along each axis for the periodic gravity mesh.
-  eta:                            0.025     # Constant dimensionless multiplier for time integration.
-  mulitpole_acceptance_criterion: adaptive  # Choice of mulitpole acceptance criterion: 'adaptive' OR 'geometric'.
-  epsilon_fmm:                    0.01      # Tolerance parameter for the adaptice multipole acceptance criterion.
-  theta_cr:                       0.7       # Opening angle for the purely gemoetric criterion.
-  comoving_DM_softening:          0.0026994 # Comoving Plummer-equivalent softening length for DM particles (in internal units).
-  max_physical_DM_softening:      0.0007    # Maximal Plummer-equivalent softening length in physical coordinates for DM particles (in internal units).
-  comoving_baryon_softening:      0.0026994 # Comoving Plummer-equivalent softening length for baryon particles (in internal units).
-  max_physical_baryon_softening:  0.0007    # Maximal Plummer-equivalent softening length in physical coordinates for baryon particles (in internal units).
-  softening_ratio_background:     0.04      # Fraction of the mean inter-particle separation to use as Plummer-equivalent softening for the background DM particles.
-  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).
-  dithering:                      1         # (Optional) Activate the dithering of the gravity mesh at every rebuild (this is the default value).
-  dithering_ratio:                1.0       # (Optional) Magnitude of each component of the dithering vector in units of the top-level cell sizes (this is the default value).
+  mesh_side_length:              128       # Number of cells along each axis for the periodic gravity mesh.
+  eta:                           0.025     # Constant dimensionless multiplier for time integration.
+  MAC:                           adaptive  # Choice of mulitpole acceptance criterion: 'adaptive' OR 'geometric'.
+  epsilon_fmm:                   0.01      # Tolerance parameter for the adaptice multipole acceptance criterion.
+  theta_cr:                      0.7       # Opening angle for the purely gemoetric criterion.
+  comoving_DM_softening:         0.0026994 # Comoving Plummer-equivalent softening length for DM particles (in internal units).
+  max_physical_DM_softening:     0.0007    # Maximal Plummer-equivalent softening length in physical coordinates for DM particles (in internal units).
+  comoving_baryon_softening:     0.0026994 # Comoving Plummer-equivalent softening length for baryon particles (in internal units).
+  max_physical_baryon_softening: 0.0007    # Maximal Plummer-equivalent softening length in physical coordinates for baryon particles (in internal units).
+  softening_ratio_background:    0.04      # Fraction of the mean inter-particle separation to use as Plummer-equivalent softening for the background DM particles.
+  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).
+  dithering:                     1         # (Optional) Activate the dithering of the gravity mesh at every rebuild (this is the default value).
+  dithering_ratio:               1.0       # (Optional) Magnitude of each component of the dithering vector in units of the top-level cell sizes (this is the default value).
 
 # Parameters when running with SWIFT_GRAVITY_FORCE_CHECKS 
 ForceChecks:
diff --git a/src/engine.c b/src/engine.c
index 1dc486e2da..2f821a95ec 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2041,19 +2041,18 @@ void engine_skip_force_and_kick(struct engine *e) {
         t->type == task_type_timestep ||
         t->type == task_type_timestep_limiter ||
         t->type == task_type_timestep_sync ||
-        t->subtype == task_subtype_force ||
-        t->subtype == task_subtype_limiter || t->subtype == task_subtype_grav ||
-        t->type == task_type_end_hydro_force ||
-        t->type == task_type_end_grav_force ||
-        t->type == task_type_grav_long_range || t->type == task_type_grav_mm ||
-        t->type == task_type_grav_down || t->type == task_type_grav_down_in ||
-        t->type == task_type_drift_gpart_out || t->type == task_type_cooling ||
+        t->type == task_type_end_hydro_force || t->type == task_type_cooling ||
         t->type == task_type_stars_in || t->type == task_type_stars_out ||
         t->type == task_type_star_formation ||
         t->type == task_type_stars_resort || t->type == task_type_extra_ghost ||
+        t->type == task_type_stars_ghost ||
+        t->type == task_type_stars_ghost_in ||
+        t->type == task_type_stars_ghost_out ||
         t->type == task_type_bh_swallow_ghost1 ||
         t->type == task_type_bh_swallow_ghost2 ||
-        t->type == task_type_bh_swallow_ghost3 ||
+        t->type == task_type_bh_swallow_ghost3 || t->type == task_type_bh_in ||
+        t->type == task_type_bh_out || t->subtype == task_subtype_force ||
+        t->subtype == task_subtype_limiter ||
         t->subtype == task_subtype_gradient ||
         t->subtype == task_subtype_stars_feedback ||
         t->subtype == task_subtype_bh_feedback ||
@@ -2069,8 +2068,7 @@ void engine_skip_force_and_kick(struct engine *e) {
         t->subtype == task_subtype_tend_gpart ||
         t->subtype == task_subtype_tend_spart ||
         t->subtype == task_subtype_tend_bpart ||
-        t->subtype == task_subtype_rho || t->subtype == task_subtype_gpart ||
-        t->subtype == task_subtype_sf_counts)
+        t->subtype == task_subtype_rho || t->subtype == task_subtype_sf_counts)
       t->skip = 1;
   }
 
@@ -2193,7 +2191,8 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
   if (e->nodeID == 0) message("Setting particles to a valid state...");
   engine_first_init_particles(e);
 
-  if (e->nodeID == 0) message("Computing initial gas densities.");
+  if (e->nodeID == 0)
+    message("Computing initial gas densities and approximate gravity.");
 
   /* Construct all cells and tasks to start everything */
   engine_rebuild(e, 0, clean_h_values);
@@ -2277,6 +2276,10 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
   space_init_sparts(e->s, e->verbose);
   space_init_bparts(e->s, e->verbose);
 
+  /* Update the MAC strategy if necessary */
+  if (e->policy & engine_policy_self_gravity)
+    gravity_props_update_MAC_choice(e->gravity_properties);
+
   /* Print the number of active tasks ? */
   if (e->verbose) engine_print_task_counts(e);
 
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index d8ce67b20b..07cfb256c3 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -87,8 +87,7 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
 
   /* Read the choice of multipole acceptance criterion */
   char buffer[32] = {0};
-  parser_get_param_string(params, "Gravity:mulitpole_acceptance_criterion",
-                          buffer);
+  parser_get_param_string(params, "Gravity:MAC", buffer);
 
   if (strcmp(buffer, "adaptive") == 0) {
     p->use_adaptive_tolerance = 1;
@@ -105,7 +104,7 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
   p->use_advanced_mac = 0;
 
   /* Geometric opening angle */
-  p->theta_crit = parser_get_param_double(params, "Gravity:theta");
+  p->theta_crit = parser_get_param_double(params, "Gravity:theta_cr");
   if (p->theta_crit >= 1.) error("Theta too large. FMM won't converge.");
 
   /* Adaptive opening angle tolerance */
@@ -190,11 +189,15 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
   gravity_props_update(p, cosmo);
 }
 
-void gravity_props_update(struct gravity_props *p,
-                          const struct cosmology *cosmo) {
+void gravity_props_update_MAC_choice(struct gravity_props *p) {
 
-  /* Choice of MAC */
+  /* Now that we have run initial accelerations,
+   * switch to the better MAC */
   if (p->use_adaptive_tolerance) p->use_advanced_mac = 1;
+}
+
+void gravity_props_update(struct gravity_props *p,
+                          const struct cosmology *cosmo) {
 
   /* Current softening length for the high-res. DM particles. */
   double DM_softening, baryon_softening;
diff --git a/src/gravity_properties.h b/src/gravity_properties.h
index 07ec862839..4f0c5433b6 100644
--- a/src/gravity_properties.h
+++ b/src/gravity_properties.h
@@ -130,7 +130,7 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
                         const int is_zoom_simulation, const int periodic);
 void gravity_props_update(struct gravity_props *p,
                           const struct cosmology *cosmo);
-
+void gravity_props_update_MAC_choice(struct gravity_props *p);
 #if defined(HAVE_HDF5)
 void gravity_props_print_snapshot(hid_t h_grpsph,
                                   const struct gravity_props *p);
-- 
GitLab