diff --git a/src/Makefile.am b/src/Makefile.am index f50f5574d16227cc925866596ad1766ee4f5a171..44ea452dc87608db483b4600607bbad500bcc51b 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -47,7 +47,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \ # Include files for distribution, not installation. nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \ - vector.h runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h \ + vector.h runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h kick.h \ gravity.h gravity_io.h \ gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \ gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \ diff --git a/src/kick.h b/src/kick.h new file mode 100644 index 0000000000000000000000000000000000000000..20ccca1eebe5a7b82ef452ae65cd00f0da539363 --- /dev/null +++ b/src/kick.h @@ -0,0 +1,110 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_KICK_H +#define SWIFT_KICK_H + +/* Config parameters. */ +#include "../config.h" + +/** + * @brief Perform the 'kick' operation on a #gpart + * + * @param gp The #gpart to kick. + * @param new_dti The (integer) time-step for this kick. + * @param timeBase The minimal allowed time-step size. + */ +__attribute__((always_inline)) INLINE static void kick_gpart(struct gpart* gp, + int new_dti, + double timeBase) { + + /* Compute the time step for this kick */ + const int ti_start = (gp->ti_begin + gp->ti_end) / 2; + const int ti_end = gp->ti_end + new_dti / 2; + const double dt = (ti_end - ti_start) * timeBase; + const double half_dt = (ti_end - gp->ti_end) * timeBase; + + /* Move particle forward in time */ + gp->ti_begin = gp->ti_end; + gp->ti_end = gp->ti_begin + new_dti; + + /* Kick particles in momentum space */ + gp->v_full[0] += gp->a_grav[0] * dt; + gp->v_full[1] += gp->a_grav[1] * dt; + gp->v_full[2] += gp->a_grav[2] * dt; + + /* Extra kick work */ + gravity_kick_extra(gp, dt, half_dt); +} + +/** + * @brief Perform the 'kick' operation on a #part + * + * @param p The #part to kick. + * @param xp The #xpart of the particle. + * @param new_dti The (integer) time-step for this kick. + * @param timeBase The minimal allowed time-step size. + */ +__attribute__((always_inline)) INLINE static void kick_part(struct part* p, + struct xpart* xp, + int new_dti, + double timeBase) { + + /* Compute the time step for this kick */ + const int ti_start = (p->ti_begin + p->ti_end) / 2; + const int ti_end = p->ti_end + new_dti / 2; + const double dt = (ti_end - ti_start) * timeBase; + const double half_dt = (ti_end - p->ti_end) * timeBase; + + /* Move particle forward in time */ + p->ti_begin = p->ti_end; + p->ti_end = p->ti_begin + new_dti; + if (p->gpart != NULL) { + p->gpart->ti_begin = p->ti_begin; + p->gpart->ti_end = p->ti_end; + } + + /* Get the acceleration */ + float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]}; + if (p->gpart != NULL) { + a_tot[0] += p->gpart->a_grav[0]; + a_tot[1] += p->gpart->a_grav[1]; + a_tot[1] += p->gpart->a_grav[2]; + } + + /* Kick particles in momentum space */ + xp->v_full[0] += a_tot[0] * dt; + xp->v_full[1] += a_tot[1] * dt; + xp->v_full[2] += a_tot[2] * dt; + if (p->gpart != NULL) { + p->gpart->v_full[0] = xp->v_full[0]; + p->gpart->v_full[1] = xp->v_full[1]; + p->gpart->v_full[2] = xp->v_full[2]; + } + + /* Go back by half-step for the hydro velocity */ + p->v[0] = xp->v_full[0] - half_dt * a_tot[0]; + p->v[1] = xp->v_full[1] - half_dt * a_tot[1]; + p->v[2] = xp->v_full[2] - half_dt * a_tot[2]; + + /* Extra kick work */ + hydro_kick_extra(p, xp, dt, half_dt); + if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt, half_dt); +} + +#endif /* SWIFT_KICK_H */ diff --git a/src/runner.c b/src/runner.c index 341370d8cb05b4f5b2721650bf59d575856c84aa..0b630b6f99bf4087ae2a7f3fd4141eddaca9a22d 100644 --- a/src/runner.c +++ b/src/runner.c @@ -47,6 +47,7 @@ #include "gravity.h" #include "hydro_properties.h" #include "hydro.h" +#include "kick.h" #include "minmax.h" #include "scheduler.h" #include "space.h" @@ -941,35 +942,20 @@ void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) { /* No children? */ if (!c->split) { - /* Loop over the g-particles and kick the active ones. */ + /* Loop over the g-particles and kick everyone. */ for (int k = 0; k < gcount; k++) { /* Get a handle on the part. */ struct gpart *const gp = &gparts[k]; - /* If the g-particle has no counterpart and needs to be kicked */ + /* If the g-particle has no counterpart */ if (gp->id < 0) { /* First, finish the force calculation */ gravity_end_force(gp); - /* Compute the time step for this kick */ - const int ti_start = (gp->ti_begin + gp->ti_end) / 2; - const int ti_end = gp->ti_end + new_dti / 2; - const double dt = (ti_end - ti_start) * timeBase; - const double half_dt = (ti_end - gp->ti_end) * timeBase; - - /* Move particle forward in time */ - gp->ti_begin = gp->ti_end; - gp->ti_end = gp->ti_begin + new_dti; - - /* Kick particles in momentum space */ - gp->v_full[0] += gp->a_grav[0] * dt; - gp->v_full[1] += gp->a_grav[1] * dt; - gp->v_full[2] += gp->a_grav[2] * dt; - - /* Extra kick work */ - gravity_kick_extra(gp, dt, half_dt); + /* Kick the g-particle forward */ + kick_gpart(gp, new_dti, timeBase); /* Number of updated g-particles */ g_updated++; @@ -982,7 +968,7 @@ void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) { /* Now do the hydro ones... */ - /* Loop over the particles and kick the active ones. */ + /* Loop over the particles and kick everyone. */ for (int k = 0; k < count; k++) { /* Get a handle on the part. */ @@ -996,47 +982,8 @@ void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) { hydro_end_force(p); if (p->gpart != NULL) gravity_end_force(p->gpart); - /* Compute the time step for this kick */ - const int ti_start = (p->ti_begin + p->ti_end) / 2; - const int ti_end = p->ti_end + new_dti / 2; - const double dt = (ti_end - ti_start) * timeBase; - const double half_dt = (ti_end - p->ti_end) * timeBase; - - /* Move particle forward in time */ - p->ti_begin = p->ti_end; - p->ti_end = p->ti_begin + new_dti; - if (p->gpart != NULL) { - p->gpart->ti_begin = p->ti_begin; - p->gpart->ti_end = p->ti_end; - } - - /* Get the acceleration */ - float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]}; - if (p->gpart != NULL) { - a_tot[0] += p->gpart->a_grav[0]; - a_tot[1] += p->gpart->a_grav[1]; - a_tot[1] += p->gpart->a_grav[2]; - } - - /* Kick particles in momentum space */ - xp->v_full[0] += a_tot[0] * dt; - xp->v_full[1] += a_tot[1] * dt; - xp->v_full[2] += a_tot[2] * dt; - - if (p->gpart != NULL) { - p->gpart->v_full[0] = xp->v_full[0]; - p->gpart->v_full[1] = xp->v_full[1]; - p->gpart->v_full[2] = xp->v_full[2]; - } - - /* Go back by half-step for the hydro velocity */ - p->v[0] = xp->v_full[0] - half_dt * a_tot[0]; - p->v[1] = xp->v_full[1] - half_dt * a_tot[1]; - p->v[2] = xp->v_full[2] - half_dt * a_tot[2]; - - /* Extra kick work */ - hydro_kick_extra(p, xp, dt, half_dt); - if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt, half_dt); + /* Kick the particle forward */ + kick_part(p, xp, new_dti, timeBase); /* Number of updated particles */ updated++; @@ -1209,24 +1156,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) { } /* Now we have a time step, proceed with the kick */ - - /* Compute the time step for this kick */ - const int ti_start = (gp->ti_begin + gp->ti_end) / 2; - const int ti_end = gp->ti_end + new_dti / 2; - const double dt = (ti_end - ti_start) * timeBase; - const double half_dt = (ti_end - gp->ti_end) * timeBase; - - /* Move particle forward in time */ - gp->ti_begin = gp->ti_end; - gp->ti_end = gp->ti_begin + new_dti; - - /* Kick particles in momentum space */ - gp->v_full[0] += gp->a_grav[0] * dt; - gp->v_full[1] += gp->a_grav[1] * dt; - gp->v_full[2] += gp->a_grav[2] * dt; - - /* Extra kick work */ - gravity_kick_extra(gp, dt, half_dt); + kick_gpart(gp, new_dti, timeBase); /* Number of updated g-particles */ g_updated++; @@ -1309,48 +1239,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) { } /* Now we have a time step, proceed with the kick */ - - /* Compute the time step for this kick */ - const int ti_start = (p->ti_begin + p->ti_end) / 2; - const int ti_end = p->ti_end + new_dti / 2; - const double dt = (ti_end - ti_start) * timeBase; - const double half_dt = (ti_end - p->ti_end) * timeBase; - - /* Move particle forward in time */ - p->ti_begin = p->ti_end; - p->ti_end = p->ti_begin + new_dti; - if (p->gpart != NULL) { - p->gpart->ti_begin = p->ti_begin; - p->gpart->ti_end = p->ti_end; - } - - /* Get the acceleration */ - float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]}; - if (p->gpart != NULL) { - a_tot[0] += p->gpart->a_grav[0]; - a_tot[1] += p->gpart->a_grav[1]; - a_tot[1] += p->gpart->a_grav[2]; - } - - /* Kick particles in momentum space */ - xp->v_full[0] += a_tot[0] * dt; - xp->v_full[1] += a_tot[1] * dt; - xp->v_full[2] += a_tot[2] * dt; - - if (p->gpart != NULL) { - p->gpart->v_full[0] = xp->v_full[0]; - p->gpart->v_full[1] = xp->v_full[1]; - p->gpart->v_full[2] = xp->v_full[2]; - } - - /* Go back by half-step for the hydro velocity */ - p->v[0] = xp->v_full[0] - half_dt * a_tot[0]; - p->v[1] = xp->v_full[1] - half_dt * a_tot[1]; - p->v[2] = xp->v_full[2] - half_dt * a_tot[2]; - - /* Extra kick work */ - hydro_kick_extra(p, xp, dt, half_dt); - if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt, half_dt); + kick_part(p, xp, new_dti, timeBase); /* Number of updated particles */ updated++;