diff --git a/src/Makefile.am b/src/Makefile.am index 8b8dce69877b1c561d96b6d5720fdab1d068436b..d9ed64cc074e124af686bc045fdfce7c0f1c23e2 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -60,7 +60,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \ kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h runner_doiact_fft.h \ units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \ - dimension.h equation_of_state.h active.h \ + dimension.h equation_of_state.h active.h timeline.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/active.h b/src/active.h index f0aba74df440e9538634bfe6a848d9ee9bcd1674..6d8a118fa892f3a339f90e19cf1565d2ca0360e1 100644 --- a/src/active.h +++ b/src/active.h @@ -27,6 +27,7 @@ #include "const.h" #include "engine.h" #include "part.h" +#include "timeline.h" /** * @brief Check that a cell been drifted to the current time. @@ -107,7 +108,8 @@ __attribute__((always_inline)) INLINE static int part_is_active( p->ti_end, e->ti_current); #endif - return (p->ti_end == e->ti_current); + // return (p->ti_end == e->ti_current); + return (p->time_bin < get_max_active_bin(e->ti_current)); } /** @@ -126,7 +128,8 @@ __attribute__((always_inline)) INLINE static int gpart_is_active( gp->ti_end, e->ti_current); #endif - return (gp->ti_end == e->ti_current); + // return (gp->ti_end == e->ti_current); + return (gp->time_bin < get_max_active_bin(e->ti_current)); } #endif /* SWIFT_ACTIVE_H */ diff --git a/src/engine.h b/src/engine.h index 861e321627032eb1f773992a9c123249c013daa0..54a0d9a3c8f19a54f87a18c8434f10d99389a1d2 100644 --- a/src/engine.h +++ b/src/engine.h @@ -82,9 +82,6 @@ extern const char *engine_policy_names[]; /* The rank of the engine as a global variable (for messages). */ extern int engine_rank; -/* The maximal number of timesteps in a simulation */ -#define max_nr_timesteps (1 << 28) - /* Data structure for the engine. */ struct engine { diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h index 9e0ca81edff06b8a32afb185f24a88b41dc87da7..4175040dc17c8d52c7b7978f138f2e911024dced 100644 --- a/src/gravity/Default/gravity.h +++ b/src/gravity/Default/gravity.h @@ -53,8 +53,7 @@ gravity_compute_timestep_self(const struct gpart* const gp) { __attribute__((always_inline)) INLINE static void gravity_first_init_gpart( struct gpart* gp) { - gp->ti_begin = 0; - gp->ti_end = 0; + gp->time_bin = -1; gp->epsilon = 0.; // MATTHIEU } diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h index f06e65e5b30ebcd609c0c6204de33da17b770add..cd16287db6a57d9e62ae4f981501208785894a7f 100644 --- a/src/gravity/Default/gravity_part.h +++ b/src/gravity/Default/gravity_part.h @@ -19,12 +19,13 @@ #ifndef SWIFT_DEFAULT_GRAVITY_PART_H #define SWIFT_DEFAULT_GRAVITY_PART_H -/* Some standard headers. */ -#include <stdlib.h> - /* Gravity particle. */ struct gpart { + /* Particle ID. If negative, it is the negative offset of the #part with + which this gpart is linked. */ + long long id_or_neg_offset; + /* Particle position. */ double x[3]; @@ -44,14 +45,13 @@ struct gpart { float epsilon; /* Particle time of beginning of time-step. */ - int ti_begin; + // int ti_begin; /* Particle time of end of time-step. */ - int ti_end; + // int ti_end; - /* Particle ID. If negative, it is the negative offset of the #part with - which this gpart is linked. */ - long long id_or_neg_offset; + /* Time-step length */ + timebin_t time_bin; } SWIFT_STRUCT_ALIGN; diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index 157893bc9e27806d2b97ac5f5a81d0f6fbb1c589..7a7d196ca324f91425018ed1523172979ca7df0e 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -212,8 +212,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep( __attribute__((always_inline)) INLINE static void hydro_first_init_part( struct part *restrict p, struct xpart *restrict xp) { - p->ti_begin = 0; - p->ti_end = 0; + p->time_bin = -1; xp->v_full[0] = p->v[0]; xp->v_full[1] = p->v[1]; xp->v_full[2] = p->v[2]; @@ -303,7 +302,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float abs_div_v = fabsf(p->density.div_v); /* Compute the pressure */ - const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase; + const integertime_t ti_begin = + get_integer_time_begin(ti_current, p->time_bin); + const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin); + const float half_dt = (ti_current - (ti_begin + ti_end) / 2) * timeBase; const float pressure = hydro_get_pressure(p, half_dt); /* Compute the sound speed */ @@ -383,7 +385,9 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( p->rho *= expf(w2); /* Drift the pressure */ - const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase; + const integertime_t ti_begin = get_integer_time_begin(t1, p->time_bin); + const integertime_t ti_end = get_integer_time_end(t1, p->time_bin); + const float dt_entr = (t1 - (ti_begin + ti_end) / 2) * timeBase; const float pressure = hydro_get_pressure(p, dt_entr); /* Compute the new sound speed */ diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h index 0c14da957463720dcee5349e88be91986cbb3f54..7c48f11b798bd79f9f515deaed6ccf99a07d6911 100644 --- a/src/hydro/Gadget2/hydro_part.h +++ b/src/hydro/Gadget2/hydro_part.h @@ -50,6 +50,12 @@ struct xpart { /* Data of a single particle. */ struct part { + /* Particle ID. */ + long long id; + + /* Pointer to corresponding gravity part. */ + struct gpart* gpart; + /* Particle position. */ double x[3]; @@ -66,10 +72,10 @@ struct part { float mass; /* Particle time of beginning of time-step. */ - int ti_begin; + // int ti_begin; /* Particle time of end of time-step. */ - int ti_end; + // int ti_end; /* Particle density. */ float rho; @@ -124,11 +130,8 @@ struct part { } force; }; - /* Particle ID. */ - long long id; - - /* Pointer to corresponding gravity part. */ - struct gpart* gpart; + /* Time-step length */ + timebin_t time_bin; } SWIFT_STRUCT_ALIGN; diff --git a/src/kick.h b/src/kick.h index e3fa3bf78c7da514abacf697a9d94212020e5a7b..91177f91b7b02b6946af03fac5c4ae7fb4a50db9 100644 --- a/src/kick.h +++ b/src/kick.h @@ -25,6 +25,7 @@ /* Local headers. */ #include "const.h" #include "debug.h" +#include "timeline.h" /** * @brief Perform the 'kick' operation on a #gpart @@ -68,8 +69,10 @@ __attribute__((always_inline)) INLINE static void kick_part( 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 integertime_t ti_begin = get_integer_time_begin(t1, p->time_bin); + const integertime_t ti_end = get_integer_time_end(t1, p->time_bin); + const int ti_start = (ti_begin + ti_end) / 2; + const int ti_end = ti_end + new_dti / 2; const float dt = (ti_end - ti_start) * timeBase; const float half_dt = (ti_end - p->ti_end) * timeBase; diff --git a/src/part.h b/src/part.h index 6832e0c6c2e0f2c324d90629447cf6a5e809d6fb..b312f1656bba8bd513d518aa5ec69ed14afc7c97 100644 --- a/src/part.h +++ b/src/part.h @@ -33,6 +33,7 @@ /* Local headers. */ #include "align.h" #include "const.h" +#include "timeline.h" /* Some constants. */ #define part_align 128 diff --git a/src/runner.c b/src/runner.c index 2d6da4e4aedc9c40d1dade243e605e9aeda86dbe..23e8aa0f698987645767dd335e69fb5ee89b5c3a 100644 --- a/src/runner.c +++ b/src/runner.c @@ -1029,7 +1029,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) { void runner_do_recv_cell(struct runner *r, struct cell *c, int timer) { const struct part *restrict parts = c->parts; - const struct gpart *restrict gparts = c->gparts; + // const struct gpart *restrict gparts = c->gparts; const size_t nr_parts = c->count; const size_t nr_gparts = c->gcount; // const int ti_current = r->e->ti_current; @@ -1045,14 +1045,14 @@ void runner_do_recv_cell(struct runner *r, struct cell *c, int timer) { /* Collect everything... */ for (size_t k = 0; k < nr_parts; k++) { - const int ti_end = parts[k].ti_end; + const int ti_end = 0; // parts[k].ti_end; //MATTHIEU // if(ti_end < ti_current) error("Received invalid particle !"); ti_end_min = min(ti_end_min, ti_end); ti_end_max = max(ti_end_max, ti_end); h_max = max(h_max, parts[k].h); } for (size_t k = 0; k < nr_gparts; k++) { - const int ti_end = gparts[k].ti_end; + const int ti_end = 0; // gparts[k].ti_end; //MATTHIEU // if(ti_end < ti_current) error("Received invalid particle !"); ti_end_min = min(ti_end_min, ti_end); ti_end_max = max(ti_end_max, ti_end); diff --git a/src/space.c b/src/space.c index e04aad2207e736ac8fe4da28b7d22b7c7999765b..0cc3c8d772199d144c95eed5630795be70ac74e8 100644 --- a/src/space.c +++ b/src/space.c @@ -1529,7 +1529,7 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data) { struct part *p = &parts[k]; struct xpart *xp = &xparts[k]; const float h = p->h; - const int ti_end = p->ti_end; + const int ti_end = get_integer_time_end(e->ti_current, p->time_bin); xp->x_diff[0] = 0.f; xp->x_diff[1] = 0.f; xp->x_diff[2] = 0.f; @@ -1539,7 +1539,7 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data) { } for (int k = 0; k < gcount; k++) { struct gpart *gp = &gparts[k]; - const int ti_end = gp->ti_end; + const int ti_end = get_integer_time_end(e->ti_current, gp->time_bin); gp->x_diff[0] = 0.f; gp->x_diff[1] = 0.f; gp->x_diff[2] = 0.f; diff --git a/src/timeline.h b/src/timeline.h new file mode 100644 index 0000000000000000000000000000000000000000..2d217a382264d178d1c969c58a67ff3c9db6c5c2 --- /dev/null +++ b/src/timeline.h @@ -0,0 +1,102 @@ +/******************************************************************************* + * 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_TIMELINE_H +#define SWIFT_TIMELINE_H + +/* Config parameters. */ +#include "../config.h" + +/* Local headers. */ +#include "inline.h" + +typedef long long integertime_t; +typedef char timebin_t; + +/*! The number of time bins */ +#define num_time_bins 20 + +/*! The maximal number of timesteps in a simulation */ +#define max_nr_timesteps (1LL << (num_time_bins + 1)) + +/** + * @brief Returns the integer time interval corresponding to a time bin + * + * @param bin The time bin of interest. + */ +static INLINE integertime_t get_integer_timestep(timebin_t bin) { + + if (bin == 0) return 0; + return 1LL << (bin + 1); +} + +/** + * @brief Returns the physical time interval corresponding to a time bin. + * + * @param bin The time bin of interest. + * @param timeBase the minimal time-step size of the simulation. + */ +static INLINE double get_timestep(timebin_t bin, double timeBase) { + + return get_integer_timestep(bin) * timeBase; +} + +/** + * @brief Returns the integer time corresponding to the start of the time-step + * given by a time-bin. + * + * @param ti_current The current time on the integer time line. + * @param bin The time bin of interest. + */ +static INLINE integertime_t get_integer_time_begin(integertime_t ti_current, + timebin_t bin) { + + const integertime_t dti = get_integer_timestep(bin); + return ti_current % dti; +} + +/** + * @brief Returns the integer time corresponding to the start of the time-step + * given by a time-bin. + * + * @param ti_current The current time on the integer time line. + * @param bin The time bin of interest. + */ +static INLINE integertime_t get_integer_time_end(integertime_t ti_current, + timebin_t bin) { + + const integertime_t dti = get_integer_timestep(bin); + return ti_current % dti + dti; +} + +/** + * @brief Returns the highest active time bin at a given point on the time line. + * + * @param time The current point on the time line. + */ +static INLINE timebin_t get_max_active_bin(integertime_t time) { + + if (time == 0) return num_time_bins; + + timebin_t bin = 1; + while (!((1LL << (bin + 1)) & time)) ++bin; + + return bin; +} + +#endif /* SWIFT_TIMELINE_H */