Recursively collect the minimal acceleration of the previous step

src/gravity/MultiSoftening/gravity.h

@@ -160,6 +160,7 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
   gp->a_grav[1] = 0.f;
   gp->a_grav[2] = 0.f;
   gp->potential = 0.f;
+  gp->old_a_grav_norm = 0.f;
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   gp->potential_PM = 0.f;
@@ -206,6 +207,14 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
   /* Apply the periodic correction to the peculiar potential */
   if (periodic) gp->potential += potential_normalisation;
 
+  /* Record the norm of the acceleration for the adaptive opening criteria.
+   * Will always be an (active) timestep behind. */
+  gp->old_a_grav_norm = gp->a_grav[0] * gp->a_grav[0] +
+                        gp->a_grav[1] * gp->a_grav[1] +
+                        gp->a_grav[2] * gp->a_grav[2];
+
+  gp->old_a_grav_norm = sqrtf(gp->old_a_grav_norm);
+
   /* Let's get physical... */
   gp->a_grav[0] *= const_G;
   gp->a_grav[1] *= const_G;

src/gravity/MultiSoftening/gravity_part.h

@@ -43,6 +43,9 @@ struct gpart {
   /*! Particle mass. */
   float mass;
 
+  /*! Norm of the acceleration at the previous step. */
+  float old_a_grav_norm;
+
   /*! Particle FoF properties (group ID, group size, ...) */
   struct fof_gpart_data fof_data;
 

src/multipole.h

@@ -52,8 +52,8 @@
 __attribute__((nonnull)) INLINE static void gravity_reset(
     struct gravity_tensors *m) {
 
-  /* Just bzero the struct. */
   bzero(m, sizeof(struct gravity_tensors));
+  m->m_pole.min_old_a_grav_norm = FLT_MAX;
 }
 
 /**
@@ -287,6 +287,7 @@ __attribute__((nonnull)) INLINE static void gravity_multipole_init(
     struct multipole *m) {
 
   bzero(m, sizeof(struct multipole));
+  m->min_old_a_grav_norm = FLT_MAX;
 }
 
 /**
@@ -356,6 +357,10 @@ __attribute__((nonnull)) INLINE static void gravity_multipole_add(
   /* Maximum of both softenings */
   ma->max_softening = max(ma->max_softening, mb->max_softening);
 
+  /* Minimum of both old accelerations */
+  ma->min_old_a_grav_norm =
+      min(ma->min_old_a_grav_norm, mb->min_old_a_grav_norm);
+
   /* Add 0th order term */
   ma->M_000 += mb->M_000;
 
@@ -482,6 +487,14 @@ __attribute__((nonnull)) INLINE static int gravity_multipole_equal(
     return 0;
   }
 
+  /* Check minimal old acceleration norm */
+  if (fabsf(ma->min_old_a_grav_norm - mb->min_old_a_grav_norm) /
+          fabsf(ma->min_old_a_grav_norm + mb->min_old_a_grav_norm) >
+      tolerance) {
+    message("min old_a_grav_norm different!");
+    return 0;
+  }
+
   /* Check bulk velocity (if non-zero and component > 1% of norm)*/
   if (fabsf(ma->vel[0] + mb->vel[0]) > 1e-10 &&
       (ma->vel[0] * ma->vel[0]) > 0.0001 * v2 &&
@@ -988,6 +1001,7 @@ __attribute__((nonnull)) INLINE static void gravity_P2M(
 
   /* Temporary variables */
   float epsilon_max = 0.f;
+  float min_old_a_grav_norm = FLT_MAX;
   double mass = 0.0;
   double com[3] = {0.0, 0.0, 0.0};
   double vel[3] = {0.f, 0.f, 0.f};
@@ -1003,6 +1017,7 @@ __attribute__((nonnull)) INLINE static void gravity_P2M(
 #endif
 
     epsilon_max = max(epsilon_max, epsilon);
+    min_old_a_grav_norm = min(min_old_a_grav_norm, gparts[k].old_a_grav_norm);
     mass += m;
     com[0] += gparts[k].x[0] * m;
     com[1] += gparts[k].x[1] * m;
@@ -1166,12 +1181,12 @@ __attribute__((nonnull)) INLINE static void gravity_P2M(
 #endif
 
   /* Store the data on the multipole. */
-  multi->m_pole.max_softening = epsilon_max;
-  multi->m_pole.M_000 = mass;
   multi->r_max = sqrt(r_max2);
   multi->CoM[0] = com[0];
   multi->CoM[1] = com[1];
   multi->CoM[2] = com[2];
+  multi->m_pole.max_softening = epsilon_max;
+  multi->m_pole.min_old_a_grav_norm = min_old_a_grav_norm;
   multi->m_pole.vel[0] = vel[0];
   multi->m_pole.vel[1] = vel[1];
   multi->m_pole.vel[2] = vel[2];
@@ -1181,6 +1196,7 @@ __attribute__((nonnull)) INLINE static void gravity_P2M(
   multi->m_pole.min_delta_vel[0] = min_delta_vel[0];
   multi->m_pole.min_delta_vel[1] = min_delta_vel[1];
   multi->m_pole.min_delta_vel[2] = min_delta_vel[2];
+  multi->m_pole.M_000 = mass;
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
 
@@ -1283,6 +1299,9 @@ __attribute__((nonnull)) INLINE static void gravity_M2M(
   /* "shift" the softening */
   m_a->max_softening = m_b->max_softening;
 
+  /* "shift" the minimal acceleration */
+  m_a->min_old_a_grav_norm = m_b->min_old_a_grav_norm;
+
   /* Shift 0th order term */
   m_a->M_000 = m_b->M_000;
 

src/multipole_struct.h

@@ -121,6 +121,9 @@ struct multipole {
   /*! Maximal co-moving softening of all the #gpart in the mulipole */
   float max_softening;
 
+  /*! Minimal acceleration norm of all the #gpart in the mulipole */
+  float min_old_a_grav_norm;
+
   /*! Mulipole power for the different orders */
   float power[SELF_GRAVITY_MULTIPOLE_ORDER + 1];