diff --git a/src/gravity/MultiSoftening/gravity.h b/src/gravity/MultiSoftening/gravity.h index d7716cf712254054e3cf364a807f3cb1c12355f3..10b956dc24e01cb1b23bf0bc6efa0c7e7e714979 100644 --- a/src/gravity/MultiSoftening/gravity.h +++ b/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; diff --git a/src/gravity/MultiSoftening/gravity_part.h b/src/gravity/MultiSoftening/gravity_part.h index 2d6d0d9bfbb18040faa282b63a0aa573b751a182..018acb2e2ddba8c0a1bbd27d7cc03a1c8633279f 100644 --- a/src/gravity/MultiSoftening/gravity_part.h +++ b/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; diff --git a/src/multipole.h b/src/multipole.h index 821c418e91b3b3ce66641ceaab3eb7a58a54249c..9827b58efe4235d296d6d6a89b914bd030cf739b 100644 --- a/src/multipole.h +++ b/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; diff --git a/src/multipole_struct.h b/src/multipole_struct.h index 907535a3239a859216bea56c010f43b8dabd1f85..f6f95576c2497aaa0d73ea109556861a90cc15cb 100644 --- a/src/multipole_struct.h +++ b/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];