From ee77845bcd2ca66bc6ec0993c07b6747a5d08f05 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Fri, 24 Jul 2020 13:58:53 +0100 Subject: [PATCH] Add fixed boundary particles to the code (boundary particles that will never have a velocity) --- configure.ac | 17 ++++ .../source/GettingStarted/special_modes.rst | 25 ++++++ src/cell.c | 2 +- src/drift.h | 83 ++++++++++++++++++- 4 files changed, 125 insertions(+), 2 deletions(-) diff --git a/configure.ac b/configure.ac index 9a287fdee..de2cbc2cf 100644 --- a/configure.ac +++ b/configure.ac @@ -409,6 +409,22 @@ elif test "$boundary_particles" != "no"; then AC_DEFINE_UNQUOTED([SWIFT_BOUNDARY_PARTICLES], [$enableval] ,[Particles with smaller ID than this will be considered as boundaries.]) fi +# Check if we want to use fixed boundary particles. +AC_ARG_ENABLE([fixed-boundary-particles], + [AS_HELP_STRING([--enable-fixed-boundary-particles], + [Set all particles with an ID smaller than @<:@N@:>@ as fixed boundary particles (i.e. receive zero gravity + hydro forces + zero velocity), this mode enables also --enable-boundary-particles and --enable-no-gravity-below-id.] + )], + [fixed_boundary_particles="$enableval"], + [fixed_boundary_particles="no"] +) +if test "$fixed_boundary_particles" == "yes"; then + AC_MSG_ERROR(Need to specify the ID below which particles get zero forces when using --enable-boundary-particles!) +elif test "$fixed_boundary_particles" != "no"; then + AC_DEFINE_UNQUOTED([SWIFT_NO_GRAVITY_BELOW_ID], [$enableval] ,[Particles with smaller ID than this will have zero gravity forces]) + AC_DEFINE_UNQUOTED([SWIFT_BOUNDARY_PARTICLES], [$enableval] ,[Particles with smaller ID than this will be considered as boundaries.]) + AC_DEFINE_UNQUOTED([SWIFT_FIXED_BOUNDARY_PARTICLES], [$enableval] ,[Particles with smaller ID than this will be considered as boundaries.]) +fi + # Check whether we have any of the ARM v8.1 tick timers AX_ASM_ARM_PMCCNTR AX_ASM_ARM_CNTVCT @@ -2395,6 +2411,7 @@ AC_MSG_RESULT([ Gravity checks : $gravity_force_checks Custom icbrtf : $enable_custom_icbrtf Boundary particles : $boundary_particles + Fixed boundary particles : $fixed_boundary_particles Particle Logger : $with_logger Python enabled : $have_python diff --git a/doc/RTD/source/GettingStarted/special_modes.rst b/doc/RTD/source/GettingStarted/special_modes.rst index ce89d857a..f9aba5db4 100644 --- a/doc/RTD/source/GettingStarted/special_modes.rst +++ b/doc/RTD/source/GettingStarted/special_modes.rst @@ -29,6 +29,31 @@ in the sampling of the halo is reduced. Note also that this does not affect the hydrodynamic forces. This mode is purely designed for gravity-only accuracy tests. +Besides the pure gravity mode, swift also has the boundary particle mode, +this mode turns off both the gravity forces and hydro forces for all +particles. Because gas particles only receive hydro this mode only impacts +gas particles more strictly than other particles. This mode can be +activated using ``--enable-boundary-particles=N``. This will zero the +gravitational and hydrodynamic *accelerations* of all particles with ``id`` +(strictly) lower than ``N`` at every time-step. Still if particles have an +initial velocity they will keep moving at that speed. This compilation +option also activates ``--enable-no-gravity-below-id=N``. + +Typical use cases are hydrodynamical experiments that have boundaries. + +Both options discussed above only cancel accelerations and have no impact +on what can happen in any code module that directly changes the velocity of +the boundary or no gravity particles. An example of this is momentum +injection of stellar winds for example. If we additionally want to keep the +boundary particles fixed at the same position for the whole simulation we can +use the ``--enable-fixed-boundary-particles=N`` compile option, this option +explicitly sets the velocity of the boundary particles to zero every time +step on top of also zeroing the accelerations. + +A typical use cases is an idealized setup in which we have a black hole in +the centre of a galaxy that accretes gas but is not allowed to move from +the momentum recoil of the gas it swallows. + Gravity glasses ~~~~~~~~~~~~~~~ diff --git a/src/cell.c b/src/cell.c index 784d9e08a..2d4b863c3 100644 --- a/src/cell.c +++ b/src/cell.c @@ -5032,7 +5032,7 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) { if (gpart_is_inhibited(gp, e)) continue; /* Drift... */ - drift_gpart(gp, dt_drift, ti_old_gpart, ti_current, grav_props); + drift_gpart(gp, dt_drift, ti_old_gpart, ti_current, grav_props, e); #ifdef SWIFT_DEBUG_CHECKS /* Make sure the particle does not drift by more than a box length. */ diff --git a/src/drift.h b/src/drift.h index 6b0a40e87..12bd716ee 100644 --- a/src/drift.h +++ b/src/drift.h @@ -27,6 +27,7 @@ #include "const.h" #include "debug.h" #include "dimension.h" +#include "engine.h" #include "entropy_floor.h" #include "hydro.h" #include "hydro_properties.h" @@ -42,10 +43,12 @@ * @param ti_old Integer start of time-step (for debugging checks). * @param ti_current Integer end of time-step (for debugging checks). * @param grav_props The properties of the gravity scheme. + * @param e the #engine */ __attribute__((always_inline)) INLINE static void drift_gpart( struct gpart *restrict gp, double dt_drift, integertime_t ti_old, - integertime_t ti_current, const struct gravity_props *grav_props) { + integertime_t ti_current, const struct gravity_props *grav_props, + const struct engine *e) { #ifdef SWIFT_DEBUG_CHECKS if (gp->ti_drift != ti_old) @@ -58,6 +61,29 @@ __attribute__((always_inline)) INLINE static void drift_gpart( gp->ti_drift = ti_current; #endif +#ifdef SWIFT_FIXED_BOUNDARY_PARTICLES + + /* Get the ID of the gpart */ + long long id = 0; + if (gp->type == swift_type_gas) + id = e->s->parts[-gp->id_or_neg_offset].id; + else if (gp->type == swift_type_stars) + id = e->s->sparts[-gp->id_or_neg_offset].id; + else if (gp->type == swift_type_black_hole) + id = e->s->bparts[-gp->id_or_neg_offset].id; + else + id = gp->id_or_neg_offset; + + /* Cancel the velocity of the particles */ + if (id < SWIFT_FIXED_BOUNDARY_PARTICLES) { + + /* Don't move! */ + gp->v_full[0] = 0.f; + gp->v_full[1] = 0.f; + gp->v_full[2] = 0.f; + } +#endif + /* Drift... */ gp->x[0] += gp->v_full[0] * dt_drift; gp->x[1] += gp->v_full[1] * dt_drift; @@ -98,6 +124,21 @@ __attribute__((always_inline)) INLINE static void drift_part( p->ti_drift = ti_current; #endif +#ifdef SWIFT_FIXED_BOUNDARY_PARTICLES + + /* Get the ID of the gpart */ + const long long id = p->id; + + /* Cancel the velocity of the particles */ + if (id < SWIFT_FIXED_BOUNDARY_PARTICLES) { + + /* Don't move! */ + xp->v_full[0] = 0.f; + xp->v_full[1] = 0.f; + xp->v_full[2] = 0.f; + } +#endif + /* Drift... */ p->x[0] += xp->v_full[0] * dt_drift; p->x[1] += xp->v_full[1] * dt_drift; @@ -121,6 +162,16 @@ __attribute__((always_inline)) INLINE static void drift_part( xp->x_diff[k] -= dx; xp->x_diff_sort[k] -= dx; } + +#ifdef SWIFT_FIXED_BOUNDARY_PARTICLES + + /* Cancel the velocity of the particles */ + if (id < SWIFT_FIXED_BOUNDARY_PARTICLES) { + p->v[0] = 0.f; + p->v[1] = 0.f; + p->v[2] = 0.f; + } +#endif } /** @@ -146,6 +197,21 @@ __attribute__((always_inline)) INLINE static void drift_spart( sp->ti_drift = ti_current; #endif +#ifdef SWIFT_FIXED_BOUNDARY_PARTICLES + + /* Get the ID of the gpart */ + const long long id = sp->id; + + /* Cancel the velocity of the particles */ + if (id < SWIFT_FIXED_BOUNDARY_PARTICLES) { + + /* Don't move! */ + sp->v[0] = 0.f; + sp->v[1] = 0.f; + sp->v[2] = 0.f; + } +#endif + /* Drift... */ sp->x[0] += sp->v[0] * dt_drift; sp->x[1] += sp->v[1] * dt_drift; @@ -185,6 +251,21 @@ __attribute__((always_inline)) INLINE static void drift_bpart( bp->ti_drift = ti_current; #endif +#ifdef SWIFT_FIXED_BOUNDARY_PARTICLES + + /* Get the ID of the gpart */ + const long long id = bp->id; + + /* Cancel the velocity of the particles */ + if (id < SWIFT_FIXED_BOUNDARY_PARTICLES) { + + /* Don't move! */ + bp->v[0] = 0.f; + bp->v[1] = 0.f; + bp->v[2] = 0.f; + } +#endif + /* Drift... */ bp->x[0] += bp->v[0] * dt_drift; bp->x[1] += bp->v[1] * dt_drift; -- 2.26.2