Commit 6d70491d authored by Matthieu Schaller's avatar Matthieu Schaller

Merge branch 'master' into test125cells

parents 97254ae5 0058b193
......@@ -31,8 +31,6 @@ Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
epsilon: 0.0001 # Softening length (in internal units).
a_smooth: 1000.
r_cut: 4.
# Parameters for the hydrodynamics scheme
SPH:
......
......@@ -26,7 +26,7 @@ from numpy import *
# with a density of 1
# Parameters
periodic= 0 # 1 For periodic box
periodic= 1 # 1 For periodic box
boxSize = 1.
rho = 1.
L = int(sys.argv[1]) # Number of particles along one axis
......
......@@ -35,4 +35,4 @@ Statistics:
# Parameters related to the initial conditions
InitialConditions:
file_name: ./uniformDMBox_100.hdf5 # The file to read
file_name: ./uniformDMBox_50.hdf5 # The file to read
......@@ -107,7 +107,7 @@ case $host_cpu in
*2?6[[ad]]?:*:*:*) ax_gcc_arch="sandybridge corei7-avx corei7 core2 pentium-m pentium3 pentiumpro" ;;
*3?6[[ae]]?:*:*:*) ax_gcc_arch="ivybridge core-avx-i corei7-avx corei7 core2 pentium-m pentium3 pentiumpro" ;;
*3?6[[cf]]?:*:*:*|*4?6[[56]]?:*:*:*) ax_gcc_arch="haswell core-avx2 core-avx-i corei7-avx corei7 core2 pentium-m pentium3 pentiumpro" ;;
*3?6d?:*:*:*) ax_gcc_arch="broadwell core-avx2 core-avx-i corei7-avx corei7 core2 pentium-m pentium3 pentiumpro" ;;
*3?6d?:*:*:*|*4?6f?:*:*:*) ax_gcc_arch="broadwell core-avx2 core-avx-i corei7-avx corei7 core2 pentium-m pentium3 pentiumpro" ;;
*1?6c?:*:*:*|*2?6[[67]]?:*:*:*|*3?6[[56]]?:*:*:*) ax_gcc_arch="bonnell atom core2 pentium-m pentium3 pentiumpro" ;;
*3?67?:*:*:*|*[[45]]?6[[ad]]?:*:*:*) ax_gcc_arch="silvermont atom core2 pentium-m pentium3 pentiumpro" ;;
*000?f[[012]]?:*:*:*|?f[[012]]?:*:*:*|f[[012]]?:*:*:*) ax_gcc_arch="pentium4 pentiumpro" ;;
......
......@@ -63,7 +63,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 cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h runner_doiact_fft.h \
runner_doiact_nosort.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 part_type.h \
dimension.h equation_of_state.h part_type.h periodic.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 \
......
......@@ -29,25 +29,48 @@
#include "timeline.h"
/**
* @brief Check that a cell been drifted to the current time.
* @brief Check that the #part in a #cell have been drifted to the current time.
*
* @param c The #cell.
* @param e The #engine containing information about the current time.
* @return 1 if the #cell has been drifted to the current time, 0 otherwise.
*/
__attribute__((always_inline)) INLINE static int cell_is_drifted(
__attribute__((always_inline)) INLINE static int cell_are_part_drifted(
const struct cell *c, const struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->ti_old > e->ti_current)
if (c->ti_old_part > e->ti_current)
error(
"Cell has been drifted too far forward in time! c->ti_old=%lld (t=%e) "
"and e->ti_current=%lld (t=%e)",
c->ti_old, c->ti_old * e->timeBase, e->ti_current,
c->ti_old_part, c->ti_old_part * e->timeBase, e->ti_current,
e->ti_current * e->timeBase);
#endif
return (c->ti_old == e->ti_current);
return (c->ti_old_part == e->ti_current);
}
/**
* @brief Check that the #gpart in a #cell have been drifted to the current
* time.
*
* @param c The #cell.
* @param e The #engine containing information about the current time.
* @return 1 if the #cell has been drifted to the current time, 0 otherwise.
*/
__attribute__((always_inline)) INLINE static int cell_are_gpart_drifted(
const struct cell *c, const struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->ti_old_gpart > e->ti_current)
error(
"Cell has been drifted too far forward in time! c->ti_old=%lld (t=%e) "
"and e->ti_current=%lld (t=%e)",
c->ti_old_gpart, c->ti_old_gpart * e->timeBase, e->ti_current,
e->ti_current * e->timeBase);
#endif
return (c->ti_old_gpart == e->ti_current);
}
/* Are cells / particles active for regular tasks ? */
......
This diff is collapsed.
......@@ -74,7 +74,7 @@ struct pcell {
/* Stats on this cell's particles. */
double h_max;
integertime_t ti_end_min, ti_end_max, ti_beg_max, ti_old;
integertime_t ti_end_min, ti_end_max, ti_beg_max, ti_old_part, ti_old_gpart;
/* Number of particles in this cell. */
int count, gcount, scount;
......@@ -157,8 +157,11 @@ struct cell {
/*! The extra ghost task for complex hydro schemes */
struct task *extra_ghost;
/*! The drift task */
struct task *drift;
/*! The drift task for parts */
struct task *drift_part;
/*! The drift task for gparts */
struct task *drift_gpart;
/*! The first kick task */
struct task *kick1;
......@@ -169,10 +172,10 @@ struct cell {
/*! The task to compute time-steps */
struct task *timestep;
/*! Task constructing the multipole from the particles */
struct task *grav_top_level;
/*! Task linking the FFT mesh to the rest of gravity tasks */
struct task *grav_ghost[2];
/*! Task constructing the multipole from the particles */
/*! Task computing long range non-periodic gravity interactions */
struct task *grav_long_range;
/*! Task propagating the multipole to the particles */
......@@ -233,24 +236,30 @@ struct cell {
/*! Maximum beginning of (integer) time step in this cell. */
integertime_t ti_beg_max;
/*! Last (integer) time the cell's particle was drifted forward in time. */
integertime_t ti_old;
/*! Last (integer) time the cell's sort arrays were updated. */
integertime_t ti_sort;
/*! Last (integer) time the cell's part were drifted forward in time. */
integertime_t ti_old_part;
/*! Last (integer) time the cell's gpart were drifted forward in time. */
integertime_t ti_old_gpart;
/*! Last (integer) time the cell's multipole was drifted forward in time. */
integertime_t ti_old_multipole;
/*! Minimum dimension, i.e. smallest edge of this cell (min(width)). */
float dmin;
/*! Maximum particle movement in this cell since last construction. */
float dx_max;
/*! Maximum particle movement in this cell since the last sort. */
float dx_max_sort;
/*! Maximum part movement in this cell since last construction. */
float dx_max_part;
/*! Maximum gpart movement in this cell since last construction. */
float dx_max_gpart;
/*! Nr of #part in this cell. */
int count;
......@@ -357,13 +366,15 @@ void cell_clean_links(struct cell *c, void *data);
void cell_make_multipoles(struct cell *c, integertime_t ti_current);
void cell_check_multipole(struct cell *c, void *data);
void cell_clean(struct cell *c);
void cell_check_particle_drift_point(struct cell *c, void *data);
void cell_check_part_drift_point(struct cell *c, void *data);
void cell_check_gpart_drift_point(struct cell *c, void *data);
void cell_check_multipole_drift_point(struct cell *c, void *data);
void cell_reset_task_counters(struct cell *c);
int cell_is_drift_needed(struct cell *c, const struct engine *e);
int cell_unskip_tasks(struct cell *c, struct scheduler *s);
void cell_set_super(struct cell *c, struct cell *super);
void cell_drift_particles(struct cell *c, const struct engine *e);
void cell_drift_part(struct cell *c, const struct engine *e);
void cell_drift_gpart(struct cell *c, const struct engine *e);
void cell_drift_multipole(struct cell *c, const struct engine *e);
void cell_drift_all_multipoles(struct cell *c, const struct engine *e);
void cell_check_timesteps(struct cell *c);
......
......@@ -74,7 +74,7 @@ hid_t io_hdf5_type(enum IO_DATA_TYPE type) {
case DOUBLE:
return H5T_NATIVE_DOUBLE;
case CHAR:
return H5T_C_S1;
return H5T_NATIVE_CHAR;
default:
error("Unknown type");
return 0;
......
......@@ -259,8 +259,8 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
result = 0;
}
if (c->dx_max != dx_max) {
message("%d Inconsistent dx_max: %f != %f", *depth, c->dx_max, dx_max);
if (c->dx_max_part != dx_max) {
message("%d Inconsistent dx_max: %f != %f", *depth, c->dx_max_part, dx_max);
message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
result = 0;
}
......
......@@ -39,7 +39,7 @@
* @param ti_current Integer end of time-step
*/
__attribute__((always_inline)) INLINE static void drift_gpart(
struct gpart *restrict gp, float dt, double timeBase, integertime_t ti_old,
struct gpart *restrict gp, double dt, double timeBase, integertime_t ti_old,
integertime_t ti_current) {
#ifdef SWIFT_DEBUG_CHECKS
......@@ -75,7 +75,7 @@ __attribute__((always_inline)) INLINE static void drift_gpart(
* @param ti_current Integer end of time-step
*/
__attribute__((always_inline)) INLINE static void drift_part(
struct part *restrict p, struct xpart *restrict xp, float dt,
struct part *restrict p, struct xpart *restrict xp, double dt,
double timeBase, integertime_t ti_old, integertime_t ti_current) {
#ifdef SWIFT_DEBUG_CHECKS
......@@ -119,7 +119,7 @@ __attribute__((always_inline)) INLINE static void drift_part(
* @param ti_current Integer end of time-step
*/
__attribute__((always_inline)) INLINE static void drift_spart(
struct spart *restrict sp, float dt, double timeBase, integertime_t ti_old,
struct spart *restrict sp, double dt, double timeBase, integertime_t ti_old,
integertime_t ti_current) {
#ifdef SWIFT_DEBUG_CHECKS
......
This diff is collapsed.
......@@ -275,7 +275,7 @@ gas_pressure_from_internal_energy(float density, float u) {
*/
__attribute__((always_inline)) INLINE static float
gas_internal_energy_from_pressure(float density, float pressure) {
return const_isothermal_energy;
return const_isothermal_internal_energy;
}
/**
......
......@@ -69,11 +69,9 @@ void gravity_props_print(const struct gravity_props *p) {
message("Self-gravity softening: epsilon=%.4f (Plummer equivalent: %.4f)",
p->epsilon, p->epsilon / 3.);
if (p->a_smooth != gravity_props_default_a_smooth)
message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth);
message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth);
if (p->r_cut != gravity_props_default_r_cut)
message("Self-gravity MM cut-off: r_cut=%f", p->r_cut);
message("Self-gravity MM cut-off: r_cut=%f", p->r_cut);
}
#if defined(HAVE_HDF5)
......
......@@ -43,18 +43,4 @@
_a > _b ? _a : _b; \
})
/**
* @brief Limits the value of x to be between a and b
*
* Only wraps once. If x > 2b, the returned value will be larger than b.
* Similarly for x < -b.
*/
#define box_wrap(x, a, b) \
({ \
const __typeof__(x) _x = (x); \
const __typeof__(a) _a = (a); \
const __typeof__(b) _b = (b); \
_x < _a ? (_x + _b) : ((_x > _b) ? (_x - _b) : _x); \
})
#endif /* SWIFT_MINMAX_H */
......@@ -28,7 +28,6 @@
#include <string.h>
/* Includes. */
//#include "active.h"
#include "align.h"
#include "const.h"
#include "error.h"
......@@ -37,8 +36,8 @@
#include "gravity_softened_derivatives.h"
#include "inline.h"
#include "kernel_gravity.h"
#include "minmax.h"
#include "part.h"
#include "periodic.h"
#include "vector_power.h"
#define multipole_align 128
......
......@@ -524,7 +524,7 @@ static void repart_edge_metis(int partweights, int bothweights, int nodeID,
t->type != task_type_sub_self && t->type != task_type_sub_self &&
t->type != task_type_ghost && t->type != task_type_kick1 &&
t->type != task_type_kick2 && t->type != task_type_timestep &&
t->type != task_type_drift)
t->type != task_type_drift_part && t->type != task_type_drift_gpart)
continue;
/* Get the task weight. */
......@@ -557,7 +557,7 @@ static void repart_edge_metis(int partweights, int bothweights, int nodeID,
/* Different weights for different tasks. */
if (t->type == task_type_ghost || t->type == task_type_kick1 ||
t->type == task_type_kick2 || t->type == task_type_timestep ||
t->type == task_type_drift) {
t->type == task_type_drift_part || t->type == task_type_drift_gpart) {
/* Particle updates add only to vertex weight. */
if (taskvweights) weights_v[cid] += w;
......
/*******************************************************************************
* 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_PERIODIC_H
#define SWIFT_PERIODIC_H
/* Config parameters. */
#include "../config.h"
/* Includes. */
#include "inline.h"
/**
* @brief Limits the value of x to be between a and b
*
* Only wraps once. If x > 2b, the returned value will be larger than b.
* Similarly for x < -b.
*/
#define box_wrap(x, a, b) \
({ \
const __typeof__(x) _x = (x); \
const __typeof__(a) _a = (a); \
const __typeof__(b) _b = (b); \
_x < _a ? (_x + _b) : ((_x > _b) ? (_x - _b) : _x); \
})
/**
* @brief Find the smallest distance dx along one axis within a box of size
* box_size
*
* This macro evaluates its arguments exactly once.
*
* Only wraps once. If dx > 2b, the returned value will be larger than b.
* Similarly for dx < -b.
*
*/
__attribute__((always_inline)) INLINE static double nearest(double dx,
double box_size) {
return dx > 0.5 * box_size ? (dx - box_size)
: ((dx < -0.5 * box_size) ? (dx + box_size) : dx);
}
/**
* @brief Find the smallest distance dx along one axis within a box of size
* box_size
*
* This macro evaluates its arguments exactly once.
*
* Only wraps once. If dx > 2b, the returned value will be larger than b.
* Similarly for dx < -b.
*
*/
__attribute__((always_inline)) INLINE static float nearestf(float dx,
float box_size) {
return dx > 0.5f * box_size
? (dx - box_size)
: ((dx < -0.5f * box_size) ? (dx + box_size) : dx);
}
#endif /* SWIFT_PERIODIC_H */
......@@ -53,6 +53,7 @@
#include "hydro_properties.h"
#include "kick.h"
#include "minmax.h"
#include "runner_doiact_fft.h"
#include "runner_doiact_vec.h"
#include "scheduler.h"
#include "sort_part.h"
......@@ -331,7 +332,7 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
TIMER_TIC;
/* Check that the particles have been moved to the current time */
if (!cell_is_drifted(c, r->e)) error("Sorting un-drifted cell");
if (!cell_are_part_drifted(c, r->e)) error("Sorting un-drifted cell");
#ifdef SWIFT_DEBUG_CHECKS
/* Make sure the sort flags are consistent (downward). */
......@@ -839,19 +840,35 @@ void runner_do_unskip_mapper(void *map_data, int num_elements,
}
}
/**
* @brief Drift particles in real space.
* @brief Drift all part in a cell.
*
* @param r The runner thread.
* @param c The cell.
* @param timer Are we timing this ?
*/
void runner_do_drift_part(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
cell_drift_part(c, r->e);
if (timer) TIMER_TOC(timer_drift_part);
}
/**
* @brief Drift all gpart in a cell.
*
* @param r The runner thread.
* @param c The cell.
* @param timer Are we timing this ?
*/
void runner_do_drift_particles(struct runner *r, struct cell *c, int timer) {
void runner_do_drift_gpart(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
cell_drift_particles(c, r->e);
cell_drift_gpart(c, r->e);
if (timer) TIMER_TOC(timer_drift);
if (timer) TIMER_TOC(timer_drift_gpart);
}
/**
......@@ -1522,7 +1539,7 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int clear_sorts,
/* ... and store. */
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
c->ti_old = ti_current;
c->ti_old_part = ti_current;
c->h_max = h_max;
if (timer) TIMER_TOC(timer_dorecv_part);
......@@ -1596,7 +1613,7 @@ void runner_do_recv_gpart(struct runner *r, struct cell *c, int timer) {
/* ... and store. */
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
c->ti_old = ti_current;
c->ti_old_gpart = ti_current;
if (timer) TIMER_TOC(timer_dorecv_gpart);
......@@ -1669,7 +1686,7 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int timer) {
/* ... and store. */
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
c->ti_old = ti_current;
c->ti_old_gpart = ti_current;
if (timer) TIMER_TOC(timer_dorecv_spart);
......@@ -1727,15 +1744,18 @@ void *runner_main(void *data) {
#ifdef SWIFT_DEBUG_CHECKS
t->ti_run = e->ti_current;
#ifndef WITH_MPI
if (ci == NULL && cj == NULL) {
if (t->type == task_type_grav_top_level) {
if (ci != NULL || cj != NULL)
error("Top-level gravity task associated with a cell");
} else if (ci == NULL && cj == NULL) {
error("Task not associated with cells!");
} else if (cj == NULL) { /* self */
if (!cell_is_active(ci, e) && t->type != task_type_sort &&
t->type != task_type_send && t->type != task_type_recv &&
t->type != task_type_kick1 && t->type != task_type_drift)
t->type != task_type_kick1 && t->type != task_type_drift_part &&
t->type != task_type_drift_gpart)
error(
"Task (type='%s/%s') should have been skipped ti_current=%lld "
"c->ti_end_min=%lld",
......@@ -1863,8 +1883,11 @@ void *runner_main(void *data) {
runner_do_extra_ghost(r, ci, 1);
break;
#endif
case task_type_drift:
runner_do_drift_particles(r, ci, 1);
case task_type_drift_part:
runner_do_drift_part(r, ci, 1);
break;
case task_type_drift_gpart:
runner_do_drift_gpart(r, ci, 1);
break;
case task_type_kick1:
runner_do_kick1(r, ci, 1);
......@@ -1902,14 +1925,11 @@ void *runner_main(void *data) {
}
break;
#endif
case task_type_grav_mm:
// runner_do_grav_mm(r, t->ci, 1);
break;
case task_type_grav_down:
runner_do_grav_down(r, t->ci, 1);
break;
case task_type_grav_top_level:
// runner_do_grav_top_level(r);
runner_do_grav_fft(r, 1);
break;
case task_type_grav_long_range:
runner_do_grav_long_range(r, t->ci, 1);
......
......@@ -62,13 +62,15 @@ struct runner {
void runner_do_ghost(struct runner *r, struct cell *c, int timer);
void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer);
void runner_do_sort(struct runner *r, struct cell *c, int flag, int clock);
void runner_do_drift_particles(struct runner *r, struct cell *c, int timer);
void runner_do_drift_part(struct runner *r, struct cell *c, int timer);
void runner_do_drift_gpart(struct runner *r, struct cell *c, int timer);
void runner_do_kick1(struct runner *r, struct cell *c, int timer);
void runner_do_kick2(struct runner *r, struct cell *c, int timer);
void runner_do_end_force(struct runner *r, struct cell *c, int timer);
void runner_do_init(struct runner *r, struct cell *c, int timer);
void runner_do_cooling(struct runner *r, struct cell *c, int timer);
void runner_do_grav_external(struct runner *r, struct cell *c, int timer);
void runner_do_grav_fft(struct runner *r, int timer);
void *runner_main(void *data);
void runner_do_unskip_mapper(void *map_data, int num_elements,
void *extra_data);
......
......@@ -903,7 +903,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
/* Anything to do here? */
if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
if (!cell_is_drifted(ci, e) || !cell_is_drifted(cj, e))
if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e))
error("Interacting undrifted cells.");
/* Get the sort ID. */
......@@ -1147,7 +1147,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
/* Anything to do here? */
if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
if (!cell_is_drifted(ci, e) || !cell_is_drifted(cj, e))
if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e))
error("Interacting undrifted cells.");
/* Get the shift ID. */
......@@ -1597,7 +1597,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
if (!cell_is_active(c, e)) return;
if (!cell_is_drifted(c, e)) error("Interacting undrifted cell.");
if (!cell_are_part_drifted(c, e)) error("Interacting undrifted cell.");
struct part *restrict parts = c->parts;
const int count = c->count;
......@@ -1846,7 +1846,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
if (!cell_is_active(c, e)) return;
if (!cell_is_drifted(c, e)) error("Cell is not drifted");
if (!cell_are_part_drifted(c, e)) error("Cell is not drifted");
struct part *restrict parts = c->parts;
const int count = c->count;
......@@ -2279,8 +2279,8 @@ void DOSUB_PAIR1(struct runner *r, struct cell *ci, struct cell *cj, int sid,
else if (cell_is_active(ci, e) || cell_is_active(cj, e)) {
/* Make sure both cells are drifted to the current timestep. */
if (!cell_is_drifted(ci, e)) cell_drift_particles(ci, e);
if (!cell_is_drifted(cj, e)) cell_drift_particles(cj, e);
if (!cell_are_part_drifted(ci, e)) cell_drift_part(ci, e);
if (!cell_are_part_drifted(cj, e)) cell_drift_part(cj, e);
/* Do any of the cells need to be sorted first? */
if (!(ci->sorted & (1 << sid)) ||
......@@ -2333,7 +2333,7 @@ void DOSUB_SELF1(struct runner *r, struct cell *ci, int gettimer) {
else {
/* Drift the cell to the current timestep if needed. */
if (!cell_is_drifted(ci, r->e)) cell_drift_particles(ci, r->e);
if (!cell_are_part_drifted(ci, r->e)) cell_drift_part(ci, r->e);
#if (DOSELF1 == runner_doself1_density) && defined(WITH_VECTORIZATION) && \
defined(GADGET2_SPH)
......@@ -2586,8 +2586,8 @@ void DOSUB_PAIR2(struct runner *r, struct cell *ci, struct cell *cj, int sid,
else if (cell_is_active(ci, e) || cell_is_active(cj, e)) {
/* Make sure both cells are drifted to the current timestep. */
if (!cell_is_drifted(ci, e)) cell_drift_particles(ci, e);
if (!cell_is_drifted(cj, e)) cell_drift_particles(cj, e);
if (!cell_are_part_drifted(ci, e)) cell_drift_part(ci, e);
if (!cell_are_part_drifted(cj, e)) cell_drift_part(cj, e);
/* Do any of the cells need to be sorted first? */
if (!(ci->sorted & (1 << sid)) ||
......@@ -3218,7 +3218,7 @@ void DOSUB_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
new_sid = sortlistID[new_sid];
/* Do any of the cells need to be drifted first? */
if (!cell_is_drifted(cj, e)) cell_drift_particles(cj, e);
if (!cell_are_part_drifted(cj, e)) cell_drift_part(cj, e);
DOPAIR_SUBSET(r, ci, parts, ind, count, cj);
}
......
......@@ -31,6 +31,283 @@
#include "runner_doiact_fft.h"
/* Local includes. */
#include "engine.h"
#include "error.h"
#include "runner.h"
#include "space.h"
#include "timers.h"
void runner_do_grav_fft(struct runner *r) {}
/**
* @brief Returns 1D index of a 3D NxNxN array using row-major style.
*
* @param i Index along x.
* @param j Index along y.
* @param k Index along z.
* @param N Size of the array along one axis.
*/
__attribute__((always_inline)) INLINE static int row_major_id(int i, int j,
int k, int N) {
return ((i % N) * N * N + (j % N) * N + (k % N));
}
/**
* @brief Assigns a given multipole to a density mesh using the CIC method.
*
* @param m The #multipole.
* @param rho The density mesh.
* @param N the size of the mesh along one axis.
* @param fac The width of a mesh cell.
*/
__attribute__((always_inline)) INLINE static void multipole_to_mesh_CIC(
const struct gravity_tensors* m, double* rho, int N, double fac) {
int i = (int)(fac * m->CoM[0]);
if (i >= N) i = N - 1;
const double dx = fac * m->CoM[0] - i;
const double tx = 1. - dx;
int j = (int)(fac * m->CoM[1]);
if (j >= N) j = N - 1;
const double dy = fac * m->CoM[1] - j;
const double ty = 1. - dy;
int k = (int)(fac * m->CoM[2]);
if (k >= N) k = N - 1;
const double dz = fac * m->CoM[2] - k;
const double tz = 1. - dz;
#ifdef SWIFT_DEBUG_CHECKS
if (i < 0 || i >= N) error("Invalid multipole position in x");
if (j < 0 || j >= N) error("Invalid multipole position in y");
if (k < 0 || k >= N) error("Invalid multipole position in z");
#endif
/* CIC ! */
rho[row_major_id(i + 0, j + 0, k + 0, N)] += m->m_pole.M_000 * tx * ty * tz;
rho[row_major_id(i + 0, j + 0, k + 1, N)] += m->m_pole.M_000 * tx * ty * dz;
rho[row_major_id(i + 0, j + 1, k + 0, N)] += m->m_pole.M_000 * tx * dy * tz;
rho[row_major_id(i + 0, j + 1, k + 1, N)] += m->m_pole.M_000 * tx * dy * dz;
rho[row_major_id(i + 1, j + 0, k + 0, N)] += m->m_pole.M_000 * dx * ty * tz;
rho[row_major_id(i + 1, j + 0, k + 1, N)] += m->m_pole.M_000 * dx * ty * dz;
rho[row_major_id(i + 1, j + 1, k + 0, N)] += m->m_pole.M_000 * dx * dy * tz;
rho[row_major_id(i + 1, j + 1, k + 1, N)] += m->m_pole.M_000 * dx * dy * dz;
}
/**
* @brief Computes the potential on a multipole from a given mesh using the CIC
* method.
*
* @param m The #multipole.
* @param pot The potential mesh.
* @param N the size of the mesh along one axis.
* @param fac width of a mesh cell.
*/
__attribute__((always_inline)) INLINE static void mesh_to_multipole_CIC(
struct gravity_tensors* m, double* pot, int N, double fac) {
int i = (int)(fac * m->CoM[0]);
if (i >= N) i = N - 1;
const double dx = fac * m->CoM[0] - i;
const double tx = 1. - dx;
int j = (int)(fac * m->CoM[1]);
if (j >= N) j = N - 1;
const double dy = fac * m->CoM[1] - j;