Skip to content
Snippets Groups Projects
Commit 3c9a3d8a authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'formatting' into 'master' 

Formatting

Fixes issue #42.

See merge request !31

Former-commit-id: e2a519b75c95e0fcf60272765868b67509cfce77
parents 4cdc697d e2bc587a
Branches
Tags
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 4161 additions and 3955 deletions
examples/test.c
+ 3
− 3
View file @ 3c9a3d8a
......@@ -795,7 +795,7 @@ int main ( int argc , char *argv[] ) {
/* Initialize the engine with this space. */
tic = getticks();
message( "nr_nodes is %i." , nr_nodes );
engine_init( &e , &s , dt_max , nr_threads , nr_queues , nr_nodes , myrank , ENGINE_POLICY | engine_policy_steal );
engine_init( &e , &s , dt_max , nr_threads , nr_queues , nr_nodes , myrank , ENGINE_POLICY | engine_policy_steal | engine_policy_paranoid );
if ( myrank == 0 )
message( "engine_init took %.3f ms." , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);
......@@ -849,12 +849,12 @@ int main ( int argc , char *argv[] ) {
/* Repartition the space amongst the nodes? */
#if defined(WITH_MPI) && defined(HAVE_METIS)
if ( j == 2 )
if ( j % 100 == 2 )
e.forcerepart = 1;
#endif
/* Force a rebuild for testing. */
/* if ( j % 4 == 1 )
/* if ( j % 4 == 3 )
e.forcerebuild = 1; */
// message( "starting run %i/%i (t=%.3e) with %i threads and %i queues..." , j+1 , runs , e.time , e.nr_threads , e.nr_queues ); fflush(stdout);
......
src/atomic.h
+ 14
− 10
View file @ 3c9a3d8a
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@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_ATOMIC_H
#define SWIFT_ATOMIC_H
/* Includes. */
#include "inline.h"
#define atomic_add(v,i) __sync_fetch_and_add( v , i )
#define atomic_inc(v) atomic_add( v , 1 )
#define atomic_dec(v) atomic_add( v , -1 )
#define atomic_cas(v,o,n) __sync_val_compare_and_swap( v , o , n )
#define atomic_add(v, i) __sync_fetch_and_add(v, i)
#define atomic_inc(v) atomic_add(v, 1)
#define atomic_dec(v) atomic_add(v, -1)
#define atomic_cas(v, o, n) __sync_val_compare_and_swap(v, o, n)
#endif /* SWIFT_ATOMIC_H */
src/cell.c
+ 445
− 470
View file @ 3c9a3d8a
This diff is collapsed.
src/cell.h
+ 141
− 133
View file @ 3c9a3d8a
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@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_CELL_H
#define SWIFT_CELL_H
/* Some constants. */
#define cell_sid_dt 13
/* Includes. */
#include "lock.h"
#include "multipole.h"
#include "part.h"
/* Forward declaration of space, needed for cell_unpack. */
struct space;
/* Some constants. */
#define cell_sid_dt 13
#define cell_max_tag (1 << 16)
/* Global variables. */
extern int cell_next_tag;
/* Packed cell. */
struct pcell {
/* Stats on this cell's particles. */
double h_max, dt_min, dt_max;
/* Number of particles in this cell. */
int count;
/* tag used for MPI communication. */
int tag;
/* Stats on this cell's particles. */
double h_max, dt_min, dt_max;
/* Number of particles in this cell. */
int count;
/* Relative indices of the cell's progeny. */
int progeny[8];
};
/* tag used for MPI communication. */
int tag;
/* Relative indices of the cell's progeny. */
int progeny[8];
};
/* Structure to store the data of a single cell. */
struct cell {
/* The cell location on the grid. */
double loc[3];
/* The cell dimensions. */
double h[3];
/* Max radii in this cell. */
double h_max;
/* Minimum and maximum dt in this cell. */
double dt_min, dt_max;
/* Minimum dimension, i.e. smallest edge of this cell. */
float dmin;
/* Maximum slack allowed for particle movement. */
float slack;
/* Maximum particle movement in this cell. */
float dx_max;
/* The depth of this cell in the tree. */
int depth, split, maxdepth;
/* Nr of parts. */
int count, gcount;
/* Pointers to the particle data. */
struct part *parts;
/* Pointers to the extra particle data. */
struct xpart *xparts;
/* Pointers to the gravity particle data. */
struct gpart *gparts;
/* Pointers for the sorted indices. */
struct entry *sort, *gsort;
unsigned int sorted, gsorted;
/* Pointers to the next level of cells. */
struct cell *progeny[8];
/* Parent cell. */
struct cell *parent;
/* Super cell, i.e. the highest-level supercell that has interactions. */
struct cell *super;
/* The task computing this cell's sorts. */
struct task *sorts, *gsorts;
int sortsize, gsortsize;
/* The tasks computing this cell's density. */
struct link *density, *force, *grav;
int nr_density, nr_force, nr_grav;
/* The ghost task to link density to interactions. */
struct task *ghost, *kick1, *kick2;
/* Task receiving data. */
struct task *recv_xv, *recv_rho;
/* Tasks for gravity tree. */
struct task *grav_up, *grav_down;
/* Number of tasks that are associated with this cell. */
int nr_tasks;
/* Is the data of this cell being used in a sub-cell? */
int hold, ghold;
/* Spin lock for various uses. */
lock_type lock, glock;
/* ID of the previous owner, e.g. runner. */
int owner;
/* Momentum of particles in cell. */
float mom[3], ang[3];
/* Potential and kinetic energy of particles in this cell. */
double epot, ekin;
/* Number of particles updated in this cell. */
int updated;
/* Linking pointer for "memory management". */
struct cell *next;
/* ID of the node this cell lives on. */
int nodeID;
/* Bit mask of the proxies this cell is registered with. */
unsigned long long int sendto;
/* Pointer to this cell's packed representation. */
struct pcell *pcell;
int pcell_size;
int tag;
/* This cell's multipole. */
struct multipole multipole;
} __attribute__((aligned (64)));
/* The cell location on the grid. */
double loc[3];
/* The cell dimensions. */
double h[3];
/* Max radii in this cell. */
double h_max;
/* Minimum and maximum dt in this cell. */
double dt_min, dt_max;
/* Minimum dimension, i.e. smallest edge of this cell. */
float dmin;
/* Maximum slack allowed for particle movement. */
float slack;
/* Maximum particle movement in this cell. */
float dx_max;
/* The depth of this cell in the tree. */
int depth, split, maxdepth;
/* Nr of parts. */
int count, gcount;
/* Pointers to the particle data. */
struct part *parts;
/* Pointers to the extra particle data. */
struct xpart *xparts;
/* Pointers to the gravity particle data. */
struct gpart *gparts;
/* Pointers for the sorted indices. */
struct entry *sort, *gsort;
unsigned int sorted, gsorted;
/* Pointers to the next level of cells. */
struct cell *progeny[8];
/* Parent cell. */
struct cell *parent;
/* Super cell, i.e. the highest-level supercell that has interactions. */
struct cell *super;
/* The task computing this cell's sorts. */
struct task *sorts, *gsorts;
int sortsize, gsortsize;
/* The tasks computing this cell's density. */
struct link *density, *force, *grav;
int nr_density, nr_force, nr_grav;
/* The ghost task to link density to interactions. */
struct task *ghost, *kick1, *kick2;
/* Task receiving data. */
struct task *recv_xv, *recv_rho;
/* Tasks for gravity tree. */
struct task *grav_up, *grav_down;
/* Number of tasks that are associated with this cell. */
int nr_tasks;
/* Is the data of this cell being used in a sub-cell? */
int hold, ghold;
/* Spin lock for various uses. */
lock_type lock, glock;
/* ID of the previous owner, e.g. runner. */
int owner;
/* Momentum of particles in cell. */
float mom[3], ang[3];
/* Potential and kinetic energy of particles in this cell. */
double epot, ekin;
/* Number of particles updated in this cell. */
int updated;
/* Linking pointer for "memory management". */
struct cell *next;
/* ID of the node this cell lives on. */
int nodeID;
/* Bit mask of the proxies this cell is registered with. */
unsigned long long int sendto;
/* Pointer to this cell's packed representation. */
struct pcell *pcell;
int pcell_size;
int tag;
/* This cell's multipole. */
struct multipole multipole;
} __attribute__((aligned(64)));
/* Function prototypes. */
void cell_split ( struct cell *c );
int cell_locktree( struct cell *c );
void cell_unlocktree( struct cell *c );
int cell_glocktree( struct cell *c );
void cell_gunlocktree( struct cell *c );
int cell_pack ( struct cell *c , struct pcell *pc );
int cell_unpack ( struct pcell *pc , struct cell *c , struct space *s );
int cell_getsize ( struct cell *c );
int cell_link ( struct cell *c , struct part *parts );
void cell_split(struct cell *c);
int cell_locktree(struct cell *c);
void cell_unlocktree(struct cell *c);
int cell_glocktree(struct cell *c);
void cell_gunlocktree(struct cell *c);
int cell_pack(struct cell *c, struct pcell *pc);
int cell_unpack(struct pcell *pc, struct cell *c, struct space *s);
int cell_getsize(struct cell *c);
int cell_link(struct cell *c, struct part *parts);
#endif /* SWIFT_CELL_H */
src/common_io.c
+ 213
− 203
View file @ 3c9a3d8a
This diff is collapsed.
src/common_io.h
+ 31
− 16
View file @ 3c9a3d8a
......@@ -2,52 +2,66 @@
* This file is part of SWIFT.
* Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
* 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_COMMON_IO_H
#define SWIFT_COMMON_IO_H
/* Config parameters. */
#include "../config.h"
/* Includes. */
#include "units.h"
#if defined(HAVE_HDF5)
/**
* @brief The different types of data used in the GADGET IC files.
*
* (This is admittedly a poor substitute to C++ templates...)
*/
enum DATA_TYPE{INT, LONG, LONGLONG, UINT, ULONG, ULONGLONG, FLOAT, DOUBLE, CHAR};
enum DATA_TYPE {
INT,
LONG,
LONGLONG,
UINT,
ULONG,
ULONGLONG,
FLOAT,
DOUBLE,
CHAR
};
/**
* @brief The two sorts of data present in the GADGET IC files: compulsory to start a run or optional.
* @brief The two sorts of data present in the GADGET IC files: compulsory to
*start a run or optional.
*
*/
enum DATA_IMPORTANCE{COMPULSORY=1, OPTIONAL=0};
enum DATA_IMPORTANCE {
COMPULSORY = 1,
OPTIONAL = 0
};
hid_t hdf5Type(enum DATA_TYPE type);
size_t sizeOfType(enum DATA_TYPE type);
void readAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data);
void writeAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data, int num);
void writeAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data,
int num);
void writeAttribute_d(hid_t grp, char* name, double data);
void writeAttribute_f(hid_t grp, char* name, float data);
......@@ -58,9 +72,9 @@ void writeAttribute_s(hid_t grp, char* name, char* str);
void createXMFfile();
FILE* prepareXMFfile();
void writeXMFfooter(FILE* xmfFile);
void writeXMFheader(FILE* xmfFile, long long N, char* hdfFileName, float time);
void writeXMFline(FILE* xmfFile, char* fileName, char* name, long long N, int dim, enum DATA_TYPE type);
void writeXMFheader(FILE* xmfFile, long long N, char* hdfFileName, float time);
void writeXMFline(FILE* xmfFile, char* fileName, char* name, long long N,
int dim, enum DATA_TYPE type);
/**
* @brief Writes the current model of SPH to the file
......@@ -74,5 +88,6 @@ void writeSPHflavour(hid_t h_file);
*/
void writeUnitSystem(hid_t h_file, struct UnitSystem* us);
#endif
#endif /* SWIFT_COMMON_IO_H */
src/const.h
+ 41
− 30
View file @ 3c9a3d8a
......@@ -2,63 +2,74 @@
* This file is part of SWIFT.
* Coypright (c) 2012 Pedro Gonnet (ptcedro.gonnet@durham.ac.uk)
* 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_CONST_H
#define SWIFT_CONST_H
/* Hydrodynamical constants. */
#define const_hydro_gamma (5.0f/3.0f)
#define const_hydro_gamma (5.0f / 3.0f)
/* SPH Viscosity constants. */
#define const_viscosity_alpha 0.8f /* Used in the legacy gadget-2 SPH mode only */
#define const_viscosity_alpha_min 0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
#define const_viscosity_alpha_max 2.0f /* Values taken from (Price,2004), not used in legacy gadget mode */
#define const_viscosity_length 0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
#define const_viscosity_alpha \
0.8f /* Used in the legacy gadget-2 SPH mode only */
#define const_viscosity_alpha_min \
0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
#define const_viscosity_alpha_max \
2.0f /* Values taken from (Price,2004), not used in legacy gadget mode */
#define const_viscosity_length \
0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
/* SPH Thermal conductivity constants. */
#define const_conductivity_alpha 1.f /* Value taken from (Price,2008), not used in legacy gadget mode */
#define const_conductivity_alpha \
1.f /* Value taken from (Price,2008), not used in legacy gadget mode */
/* Time integration constants. */
#define const_cfl 0.3f
#define const_ln_max_h_change 0.231111721f /* Particle can't change volume by more than a factor of 2=1.26^3 over one time step */
#define const_max_u_change 0.1f
#define const_cfl 0.3f
#define const_ln_max_h_change \
0.231111721f /* Particle can't change volume by more than a factor of \
2=1.26^3 over one time step */
#define const_max_u_change 0.1f
/* Neighbour search constants. */
#define const_eta_kernel 1.2349f /* Corresponds to 48 ngbs with the cubic spline kernel */
#define const_delta_nwneigh 1.f
#define const_eta_kernel \
1.2349f /* Corresponds to 48 ngbs with the cubic spline kernel */
#define const_delta_nwneigh 1.f
#define CUBIC_SPLINE_KERNEL
/* Gravity stuff. */
#define const_theta_max 0.57735f /* Opening criteria, which is the ratio of the
cell distance over the cell width. */
#define const_theta_max \
0.57735f /* Opening criteria, which is the ratio of the \
cell distance over the cell width. */
// #define const_G 6.67384e-8f /* Gravitational constant. */
#define const_G 6.672e-8f /* Gravitational constant. */
#define const_epsilon 0.0014f /* Gravity blending distance. */
#define const_iepsilon 714.285714286f /* Inverse gravity blending distance. */
#define const_iepsilon2 (const_iepsilon*const_iepsilon)
#define const_iepsilon3 (const_iepsilon2*const_iepsilon)
#define const_iepsilon4 (const_iepsilon2*const_iepsilon2)
#define const_iepsilon5 (const_iepsilon3*const_iepsilon2)
#define const_iepsilon6 (const_iepsilon3*const_iepsilon3)
#define const_G 6.672e-8f /* Gravitational constant. */
#define const_epsilon 0.0014f /* Gravity blending distance. */
#define const_iepsilon 714.285714286f /* Inverse gravity blending distance. */
#define const_iepsilon2 (const_iepsilon* const_iepsilon)
#define const_iepsilon3 (const_iepsilon2* const_iepsilon)
#define const_iepsilon4 (const_iepsilon2* const_iepsilon2)
#define const_iepsilon5 (const_iepsilon3* const_iepsilon2)
#define const_iepsilon6 (const_iepsilon3* const_iepsilon3)
/* SPH variant to use */
#define LEGACY_GADGET2_SPH
/* System of units */
#define const_unit_length_in_cgs 1 /* 3.08567810e16 /\* 1Mpc *\/ */
#define const_unit_mass_in_cgs 1 /* 1.9891e33 /\* 1 M_sun *\/ */
#define const_unit_velocity_in_cgs 1 /* 1e5 /\* km s^-1 *\/ */
#define const_unit_length_in_cgs 1 /* 3.08567810e16 /\* 1Mpc *\/ */
#define const_unit_mass_in_cgs 1 /* 1.9891e33 /\* 1 M_sun *\/ */
#define const_unit_velocity_in_cgs 1 /* 1e5 /\* km s^-1 *\/ */
#endif /* SWIFT_CONST_H */
src/cycle.h
+ 156
− 165
View file @ 3c9a3d8a
......@@ -23,7 +23,6 @@
*
*/
/* machine-dependent cycle counters code. Needs to be inlined. */
/***************************************************************************/
......@@ -52,25 +51,28 @@
defined according to whether the corresponding function/type/header
is available on your system. The necessary macros are most
conveniently defined if you are using GNU autoconf, via the tests:
dnl ---------------------------------------------------------------------
AC_C_INLINE
AC_HEADER_TIME
AC_CHECK_HEADERS([sys/time.h c_asm.h intrinsics.h mach/mach_time.h])
AC_CHECK_TYPE([hrtime_t],[AC_DEFINE(HAVE_HRTIME_T, 1, [Define to 1 if hrtime_t is defined in <sys/time.h>])],,[#if HAVE_SYS_TIME_H
AC_CHECK_TYPE([hrtime_t],[AC_DEFINE(HAVE_HRTIME_T, 1, [Define to 1 if
hrtime_t is defined in <sys/time.h>])],,[#if HAVE_SYS_TIME_H
#include <sys/time.h>
#endif])
AC_CHECK_FUNCS([gethrtime read_real_time time_base_to_time clock_gettime mach_absolute_time])
AC_CHECK_FUNCS([gethrtime read_real_time time_base_to_time clock_gettime
mach_absolute_time])
dnl Cray UNICOS _rtc() (real-time clock) intrinsic
AC_MSG_CHECKING([for _rtc intrinsic])
rtc_ok=yes
AC_TRY_LINK([#ifdef HAVE_INTRINSICS_H
#include <intrinsics.h>
#endif], [_rtc()], [AC_DEFINE(HAVE__RTC,1,[Define if you have the UNICOS _rtc() intrinsic.])], [rtc_ok=no])
#endif], [_rtc()], [AC_DEFINE(HAVE__RTC,1,[Define if you have the UNICOS _rtc()
intrinsic.])], [rtc_ok=no])
AC_MSG_RESULT($rtc_ok)
dnl ---------------------------------------------------------------------
......@@ -79,24 +81,25 @@
/***************************************************************************/
#if TIME_WITH_SYS_TIME
# include <sys/time.h>
# include <time.h>
#include <sys/time.h>
#include <time.h>
#else
# if HAVE_SYS_TIME_H
# include <sys/time.h>
# else
# include <time.h>
# endif
#if HAVE_SYS_TIME_H
#include <sys/time.h>
#else
#include <time.h>
#endif
#endif
#define INLINE_ELAPSED(INL) static INL double elapsed(ticks t1, ticks t0) \
{ \
return (double)t1 - (double)t0; \
}
#define INLINE_ELAPSED(INL) \
static INL double elapsed(ticks t1, ticks t0) { \
return (double)t1 - (double)t0; \
}
/*----------------------------------------------------------------*/
/* Solaris */
#if defined(HAVE_GETHRTIME) && defined(HAVE_HRTIME_T) && !defined(HAVE_TICK_COUNTER)
#if defined(HAVE_GETHRTIME) && defined(HAVE_HRTIME_T) && \
!defined(HAVE_TICK_COUNTER)
typedef hrtime_t ticks;
#define getticks gethrtime
......@@ -108,22 +111,22 @@ INLINE_ELAPSED(inline)
/*----------------------------------------------------------------*/
/* AIX v. 4+ routines to read the real-time clock or time-base register */
#if defined(HAVE_READ_REAL_TIME) && defined(HAVE_TIME_BASE_TO_TIME) && !defined(HAVE_TICK_COUNTER)
#if defined(HAVE_READ_REAL_TIME) && defined(HAVE_TIME_BASE_TO_TIME) && \
!defined(HAVE_TICK_COUNTER)
typedef timebasestruct_t ticks;
static __inline ticks getticks(void)
{
ticks t;
read_real_time(&t, TIMEBASE_SZ);
return t;
static __inline ticks getticks(void) {
ticks t;
read_real_time(&t, TIMEBASE_SZ);
return t;
}
static __inline double elapsed(ticks t1, ticks t0) /* time in nanoseconds */
{
time_base_to_time(&t1, TIMEBASE_SZ);
time_base_to_time(&t0, TIMEBASE_SZ);
return (((double)t1.tb_high - (double)t0.tb_high) * 1.0e9 +
((double)t1.tb_low - (double)t0.tb_low));
time_base_to_time(&t1, TIMEBASE_SZ);
time_base_to_time(&t0, TIMEBASE_SZ);
return (((double)t1.tb_high - (double)t0.tb_high) * 1.0e9 +
((double)t1.tb_low - (double)t0.tb_low));
}
#define HAVE_TICK_COUNTER
......@@ -133,20 +136,23 @@ static __inline double elapsed(ticks t1, ticks t0) /* time in nanoseconds */
/*
* PowerPC ``cycle'' counter using the time base register.
*/
#if ((((defined(__GNUC__) && (defined(__powerpc__) || defined(__ppc__))) || (defined(__MWERKS__) && defined(macintosh)))) || (defined(__IBM_GCC_ASM) && (defined(__powerpc__) || defined(__ppc__)))) && !defined(HAVE_TICK_COUNTER)
#if ((((defined(__GNUC__) && (defined(__powerpc__) || defined(__ppc__))) || \
(defined(__MWERKS__) && defined(macintosh)))) || \
(defined(__IBM_GCC_ASM) && \
(defined(__powerpc__) || defined(__ppc__)))) && \
!defined(HAVE_TICK_COUNTER)
typedef unsigned long long ticks;
static __inline__ ticks getticks(void)
{
unsigned int tbl, tbu0, tbu1;
static __inline__ ticks getticks(void) {
unsigned int tbl, tbu0, tbu1;
do {
__asm__ __volatile__ ("mftbu %0" : "=r"(tbu0));
__asm__ __volatile__ ("mftb %0" : "=r"(tbl));
__asm__ __volatile__ ("mftbu %0" : "=r"(tbu1));
} while (tbu0 != tbu1);
do {
__asm__ __volatile__("mftbu %0" : "=r"(tbu0));
__asm__ __volatile__("mftb %0" : "=r"(tbl));
__asm__ __volatile__("mftbu %0" : "=r"(tbu1));
} while (tbu0 != tbu1);
return (((unsigned long long)tbu0) << 32) | tbl;
return (((unsigned long long)tbu0) << 32) | tbl;
}
INLINE_ELAPSED(__inline__)
......@@ -156,7 +162,8 @@ INLINE_ELAPSED(__inline__)
/* MacOS/Mach (Darwin) time-base register interface (unlike UpTime,
from Carbon, requires no additional libraries to be linked). */
#if defined(HAVE_MACH_ABSOLUTE_TIME) && defined(HAVE_MACH_MACH_TIME_H) && !defined(HAVE_TICK_COUNTER)
#if defined(HAVE_MACH_ABSOLUTE_TIME) && defined(HAVE_MACH_MACH_TIME_H) && \
!defined(HAVE_TICK_COUNTER)
#include <mach/mach_time.h>
typedef uint64_t ticks;
#define getticks mach_absolute_time
......@@ -166,31 +173,31 @@ INLINE_ELAPSED(__inline__)
/*----------------------------------------------------------------*/
/*
* Pentium cycle counter
* Pentium cycle counter
*/
#if (defined(__GNUC__) || defined(__ICC)) && defined(__i386__) && !defined(HAVE_TICK_COUNTER)
#if (defined(__GNUC__) || defined(__ICC)) && defined(__i386__) && \
!defined(HAVE_TICK_COUNTER)
typedef unsigned long long ticks;
#ifndef INLINE
# if __GNUC__ && !__GNUC_STDC_INLINE__
# define INLINE extern inline
# else
# define INLINE inline
# endif
#if __GNUC__ && !__GNUC_STDC_INLINE__
#define INLINE extern inline
#else
#define INLINE inline
#endif
INLINE static ticks getticks(void)
{
ticks ret;
#endif
INLINE static ticks getticks(void) {
ticks ret;
__asm__ __volatile__("rdtsc": "=A" (ret));
/* no input, nothing else clobbered */
return ret;
__asm__ __volatile__("rdtsc" : "=A"(ret));
/* no input, nothing else clobbered */
return ret;
}
INLINE_ELAPSED(__inline__)
#define HAVE_TICK_COUNTER
#define TIME_MIN 5000.0 /* unreliable pentium IV cycle counter */
#define TIME_MIN 5000.0 /* unreliable pentium IV cycle counter */
#endif
/* Visual C++ -- thanks to Morten Nissov for his help with this */
......@@ -199,46 +206,43 @@ INLINE_ELAPSED(__inline__)
typedef LARGE_INTEGER ticks;
#define RDTSC __asm __emit 0fh __asm __emit 031h /* hack for VC++ 5.0 */
static __inline ticks getticks(void)
{
ticks retval;
__asm {
RDTSC
mov retval.HighPart, edx
mov retval.LowPart, eax
}
return retval;
static __inline ticks getticks(void) {
ticks retval;
__asm {
RDTSC
mov retval.HighPart, edx mov retval.LowPart, eax
}
return retval;
}
static __inline double elapsed(ticks t1, ticks t0)
{
return (double)t1.QuadPart - (double)t0.QuadPart;
}
static __inline double elapsed(ticks t1, ticks t0) {
return (double)t1.QuadPart - (double)t0.QuadPart;
}
#define HAVE_TICK_COUNTER
#define TIME_MIN 5000.0 /* unreliable pentium IV cycle counter */
#define TIME_MIN 5000.0 /* unreliable pentium IV cycle counter */
#endif
/*----------------------------------------------------------------*/
/*
* X86-64 cycle counter
*/
#if (defined(__GNUC__) || defined(__ICC) || defined(__SUNPRO_C)) && defined(__x86_64__) && !defined(HAVE_TICK_COUNTER)
#if (defined(__GNUC__) || defined(__ICC) || defined(__SUNPRO_C)) && \
defined(__x86_64__) && !defined(HAVE_TICK_COUNTER)
typedef unsigned long long ticks;
#ifndef INLINE
# if __GNUC__ && !__GNUC_STDC_INLINE__
# define INLINE extern inline
# else
# define INLINE inline
# endif
#if __GNUC__ && !__GNUC_STDC_INLINE__
#define INLINE extern inline
#else
#define INLINE inline
#endif
INLINE static ticks getticks(void)
{
unsigned a, d;
asm volatile("rdtsc" : "=a" (a), "=d" (d));
return ((ticks)a) | (((ticks)d) << 32);
#endif
INLINE static ticks getticks(void) {
unsigned a, d;
asm volatile("rdtsc" : "=a"(a), "=d"(d));
return ((ticks)a) | (((ticks)d) << 32);
}
INLINE_ELAPSED(__inline__)
......@@ -249,18 +253,18 @@ INLINE_ELAPSED(__inline__)
/* PGI compiler, courtesy Cristiano Calonaci, Andrea Tarsi, & Roberto Gori.
NOTE: this code will fail to link unless you use the -Masmkeyword compiler
option (grrr). */
#if defined(__PGI) && defined(__x86_64__) && !defined(HAVE_TICK_COUNTER)
#if defined(__PGI) && defined(__x86_64__) && !defined(HAVE_TICK_COUNTER)
typedef unsigned long long ticks;
static ticks getticks(void)
{
asm(" rdtsc; shl $0x20,%rdx; mov %eax,%eax; or %rdx,%rax; ");
static ticks getticks(void) {
asm(" rdtsc; shl $0x20,%rdx; mov %eax,%eax; or %rdx,%rax; ");
}
INLINE_ELAPSED(__inline__)
#define HAVE_TICK_COUNTER
#endif
/* Visual C++, courtesy of Dirk Michaelis */
#if _MSC_VER >= 1400 && (defined(_M_AMD64) || defined(_M_X64)) && !defined(HAVE_TICK_COUNTER)
#if _MSC_VER >= 1400 && (defined(_M_AMD64) || defined(_M_X64)) && \
!defined(HAVE_TICK_COUNTER)
#include <intrin.h>
#pragma intrinsic(__rdtsc)
......@@ -277,17 +281,15 @@ INLINE_ELAPSED(__inline)
*/
/* intel's icc/ecc compiler */
#if (defined(__EDG_VERSION) || defined(__ECC)) && defined(__ia64__) && !defined(HAVE_TICK_COUNTER)
#if (defined(__EDG_VERSION) || defined(__ECC)) && defined(__ia64__) && \
!defined(HAVE_TICK_COUNTER)
typedef unsigned long ticks;
#include <ia64intrin.h>
static __inline__ ticks getticks(void)
{
return __getReg(_IA64_REG_AR_ITC);
}
static __inline__ ticks getticks(void) { return __getReg(_IA64_REG_AR_ITC); }
INLINE_ELAPSED(__inline__)
#define HAVE_TICK_COUNTER
#endif
......@@ -295,12 +297,11 @@ INLINE_ELAPSED(__inline__)
#if defined(__GNUC__) && defined(__ia64__) && !defined(HAVE_TICK_COUNTER)
typedef unsigned long ticks;
static __inline__ ticks getticks(void)
{
ticks ret;
static __inline__ ticks getticks(void) {
ticks ret;
__asm__ __volatile__ ("mov %0=ar.itc" : "=r"(ret));
return ret;
__asm__ __volatile__("mov %0=ar.itc" : "=r"(ret));
return ret;
}
INLINE_ELAPSED(__inline__)
......@@ -313,12 +314,11 @@ INLINE_ELAPSED(__inline__)
#include <machine/sys/inline.h>
typedef unsigned long ticks;
static inline ticks getticks(void)
{
ticks ret;
static inline ticks getticks(void) {
ticks ret;
ret = _Asm_mov_from_ar (_AREG_ITC);
return ret;
ret = _Asm_mov_from_ar(_AREG_ITC);
return ret;
}
INLINE_ELAPSED(inline)
......@@ -330,17 +330,17 @@ INLINE_ELAPSED(inline)
#if defined(_MSC_VER) && defined(_M_IA64) && !defined(HAVE_TICK_COUNTER)
typedef unsigned __int64 ticks;
# ifdef __cplusplus
#ifdef __cplusplus
extern "C"
# endif
ticks __getReg(int whichReg);
#endif
ticks
__getReg(int whichReg);
#pragma intrinsic(__getReg)
static __inline ticks getticks(void)
{
volatile ticks temp;
temp = __getReg(3116);
return temp;
static __inline ticks getticks(void) {
volatile ticks temp;
temp = __getReg(3116);
return temp;
}
INLINE_ELAPSED(inline)
......@@ -350,29 +350,27 @@ INLINE_ELAPSED(inline)
/*----------------------------------------------------------------*/
/*
* PA-RISC cycle counter
* PA-RISC cycle counter
*/
#if defined(__hppa__) || defined(__hppa) && !defined(HAVE_TICK_COUNTER)
typedef unsigned long ticks;
# ifdef __GNUC__
static __inline__ ticks getticks(void)
{
ticks ret;
#ifdef __GNUC__
static __inline__ ticks getticks(void) {
ticks ret;
__asm__ __volatile__("mfctl 16, %0": "=r" (ret));
/* no input, nothing else clobbered */
return ret;
__asm__ __volatile__("mfctl 16, %0" : "=r"(ret));
/* no input, nothing else clobbered */
return ret;
}
# else
# include <machine/inline.h>
static inline unsigned long getticks(void)
{
register ticks ret;
_MFCTL(16, ret);
return ret;
#else
#include <machine/inline.h>
static inline unsigned long getticks(void) {
register ticks ret;
_MFCTL(16, ret);
return ret;
}
# endif
#endif
INLINE_ELAPSED(inline)
......@@ -384,11 +382,10 @@ INLINE_ELAPSED(inline)
#if defined(__GNUC__) && defined(__s390__) && !defined(HAVE_TICK_COUNTER)
typedef unsigned long long ticks;
static __inline__ ticks getticks(void)
{
ticks cycles;
__asm__("stck 0(%0)" : : "a" (&(cycles)) : "memory", "cc");
return cycles;
static __inline__ ticks getticks(void) {
ticks cycles;
__asm__("stck 0(%0)" : : "a"(&(cycles)) : "memory", "cc");
return cycles;
}
INLINE_ELAPSED(__inline__)
......@@ -398,16 +395,15 @@ INLINE_ELAPSED(__inline__)
/*----------------------------------------------------------------*/
#if defined(__GNUC__) && defined(__alpha__) && !defined(HAVE_TICK_COUNTER)
/*
* The 32-bit cycle counter on alpha overflows pretty quickly,
* The 32-bit cycle counter on alpha overflows pretty quickly,
* unfortunately. A 1GHz machine overflows in 4 seconds.
*/
typedef unsigned int ticks;
static __inline__ ticks getticks(void)
{
unsigned long cc;
__asm__ __volatile__ ("rpcc %0" : "=r"(cc));
return (cc & 0xFFFFFFFF);
static __inline__ ticks getticks(void) {
unsigned long cc;
__asm__ __volatile__("rpcc %0" : "=r"(cc));
return (cc & 0xFFFFFFFF);
}
INLINE_ELAPSED(__inline__)
......@@ -419,11 +415,10 @@ INLINE_ELAPSED(__inline__)
#if defined(__GNUC__) && defined(__sparc_v9__) && !defined(HAVE_TICK_COUNTER)
typedef unsigned long ticks;
static __inline__ ticks getticks(void)
{
ticks ret;
__asm__ __volatile__("rd %%tick, %0" : "=r" (ret));
return ret;
static __inline__ ticks getticks(void) {
ticks ret;
__asm__ __volatile__("rd %%tick, %0" : "=r"(ret));
return ret;
}
INLINE_ELAPSED(__inline__)
......@@ -432,15 +427,15 @@ INLINE_ELAPSED(__inline__)
#endif
/*----------------------------------------------------------------*/
#if (defined(__DECC) || defined(__DECCXX)) && defined(__alpha) && defined(HAVE_C_ASM_H) && !defined(HAVE_TICK_COUNTER)
# include <c_asm.h>
#if (defined(__DECC) || defined(__DECCXX)) && defined(__alpha) && \
defined(HAVE_C_ASM_H) && !defined(HAVE_TICK_COUNTER)
#include <c_asm.h>
typedef unsigned int ticks;
static __inline ticks getticks(void)
{
unsigned long cc;
cc = asm("rpcc %v0");
return (cc & 0xFFFFFFFF);
static __inline ticks getticks(void) {
unsigned long cc;
cc = asm("rpcc %v0");
return (cc & 0xFFFFFFFF);
}
INLINE_ELAPSED(__inline)
......@@ -449,20 +444,19 @@ INLINE_ELAPSED(__inline)
#endif
/*----------------------------------------------------------------*/
/* SGI/Irix */
#if defined(HAVE_CLOCK_GETTIME) && defined(CLOCK_SGI_CYCLE) && !defined(HAVE_TICK_COUNTER)
#if defined(HAVE_CLOCK_GETTIME) && defined(CLOCK_SGI_CYCLE) && \
!defined(HAVE_TICK_COUNTER)
typedef struct timespec ticks;
static inline ticks getticks(void)
{
struct timespec t;
clock_gettime(CLOCK_SGI_CYCLE, &t);
return t;
static inline ticks getticks(void) {
struct timespec t;
clock_gettime(CLOCK_SGI_CYCLE, &t);
return t;
}
static inline double elapsed(ticks t1, ticks t0)
{
return ((double)t1.tv_sec - (double)t0.tv_sec) * 1.0E9 +
((double)t1.tv_nsec - (double)t0.tv_nsec);
static inline double elapsed(ticks t1, ticks t0) {
return ((double)t1.tv_sec - (double)t0.tv_sec) * 1.0E9 +
((double)t1.tv_nsec - (double)t0.tv_nsec);
}
#define HAVE_TICK_COUNTER
#endif
......@@ -471,7 +465,7 @@ static inline double elapsed(ticks t1, ticks t0)
/* Cray UNICOS _rtc() intrinsic function */
#if defined(HAVE__RTC) && !defined(HAVE_TICK_COUNTER)
#ifdef HAVE_INTRINSICS_H
# include <intrinsics.h>
#include <intrinsics.h>
#endif
typedef long long ticks;
......@@ -493,25 +487,23 @@ INLINE_ELAPSED(inline)
typedef uint64_t ticks;
static inline ticks getticks(void)
{
static inline ticks getticks(void) {
static uint64_t* addr = 0;
if (addr == 0)
{
if (addr == 0) {
uint32_t rq_addr = 0x10030000;
int fd;
int pgsize;
pgsize = getpagesize();
fd = open ("/dev/mem", O_RDONLY | O_SYNC, 0);
fd = open("/dev/mem", O_RDONLY | O_SYNC, 0);
if (fd < 0) {
perror("open");
return NULL;
}
addr = mmap(0, pgsize, PROT_READ, MAP_SHARED, fd, rq_addr);
close(fd);
if (addr == (uint64_t *)-1) {
if (addr == (uint64_t*)-1) {
perror("mmap");
return NULL;
}
......@@ -525,4 +517,3 @@ INLINE_ELAPSED(inline)
#define HAVE_TICK_COUNTER
#endif
#endif /* HAVE_MIPS_ZBUS_TIMER */
src/debug.c
+ 74
− 88
View file @ 3c9a3d8a
......@@ -2,125 +2,111 @@
* This file is part of SWIFT.
* Coypright (c) 2013 Matthieu Schaller (matthieu.schaller@durham.ac.uk),
* Pedro Gonnet (pedro.gonnet@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/>.
*
*
******************************************************************************/
#include <stdio.h>
#include "const.h"
#include "part.h"
/* This object's header. */
#include "debug.h"
/**
* @brief Dump the information pertaining to the given cell.
*/
void print_cell(struct cell *c) {
printf(
"## Cell 0x%0zx: loc=[%.3e,%.3e,%.3e], h=[%.3e,%.3e,%.3e], depth=%i, "
"split=%i, maxdepth=%i.\n",
(size_t)c, c->loc[0], c->loc[1], c->loc[2], c->h[0], c->h[1], c->h[2],
c->depth, c->split, c->maxdepth);
}
/**
* @brief Looks for the particle with the given id and prints its information to the standard output.
*
* @brief Looks for the particle with the given id and prints its information to
*the standard output.
*
* @param parts The array of particles.
* @param id The id too look for.
* @param N The size of the array of particles.
*
* (Should be used for debugging only as it runs in O(N).)
*/
void printParticle ( struct part *parts , long long int id, int N ) {
int i, found = 0;
/* Look for the particle. */
for ( i = 0 ; i < N ; i++ )
if ( parts[i].id == id ) {
printf("## Particle[%d]: id=%lld, x=[%.16e,%.16e,%.16e], v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], h=%.3e, h_dt=%.3e, wcount=%.3e, m=%.3e, rho=%.3e, rho_dh=%.3e, div_v=%.3e, u=%.3e, dudt=%.3e, bals=%.3e, POrho2=%.3e, v_sig=%.3e, dt=%.3e\n",
i,
parts[i].id,
parts[i].x[0], parts[i].x[1], parts[i].x[2],
parts[i].v[0], parts[i].v[1], parts[i].v[2],
parts[i].a[0], parts[i].a[1], parts[i].a[2],
parts[i].h,
parts[i].force.h_dt,
parts[i].density.wcount,
parts[i].mass,
parts[i].rho, parts[i].rho_dh,
parts[i].density.div_v,
parts[i].u,
parts[i].force.u_dt,
parts[i].force.balsara,
parts[i].force.POrho2,
parts[i].force.v_sig,
parts[i].dt
);
found = 1;
}
if ( !found )
printf("## Particles[???] id=%lld not found\n", id);
void printParticle(struct part *parts, long long int id, int N) {
int i, found = 0;
/* Look for the particle. */
for (i = 0; i < N; i++)
if (parts[i].id == id) {
printf(
"## Particle[%d]: id=%lld, x=[%.16e,%.16e,%.16e], "
"v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], h=%.3e, h_dt=%.3e, "
"wcount=%.3e, m=%.3e, rho=%.3e, rho_dh=%.3e, div_v=%.3e, u=%.3e, "
"dudt=%.3e, bals=%.3e, POrho2=%.3e, v_sig=%.3e, dt=%.3e\n",
i, parts[i].id, parts[i].x[0], parts[i].x[1], parts[i].x[2],
parts[i].v[0], parts[i].v[1], parts[i].v[2], parts[i].a[0],
parts[i].a[1], parts[i].a[2], parts[i].h, parts[i].force.h_dt,
parts[i].density.wcount, parts[i].mass, parts[i].rho, parts[i].rho_dh,
parts[i].density.div_v, parts[i].u, parts[i].force.u_dt,
parts[i].force.balsara, parts[i].force.POrho2, parts[i].force.v_sig,
parts[i].dt);
found = 1;
}
if (!found) printf("## Particles[???] id=%lld not found\n", id);
}
void printgParticle(struct gpart *parts, long long int id, int N) {
int i, found = 0;
void printgParticle ( struct gpart *parts , long long int id, int N ) {
int i, found = 0;
/* Look for the particle. */
for ( i = 0 ; i < N ; i++ )
if ( parts[i].id == -id || ( parts[i].id > 0 && parts[i].part->id == id ) ) {
printf("## gParticle[%d]: id=%lld, x=[%.16e,%.16e,%.16e], v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], m=%.3e, dt=%.3e\n",
i,
(parts[i].id < 0) ? -parts[i].id : parts[i].part->id ,
parts[i].x[0], parts[i].x[1], parts[i].x[2],
parts[i].v[0], parts[i].v[1], parts[i].v[2],
parts[i].a[0], parts[i].a[1], parts[i].a[2],
parts[i].mass,
parts[i].dt
);
found = 1;
}
if ( !found )
printf("## Particles[???] id=%lld not found\n", id);
/* Look for the particle. */
for (i = 0; i < N; i++)
if (parts[i].id == -id || (parts[i].id > 0 && parts[i].part->id == id)) {
printf(
"## gParticle[%d]: id=%lld, x=[%.16e,%.16e,%.16e], "
"v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], m=%.3e, dt=%.3e\n",
i, (parts[i].id < 0) ? -parts[i].id : parts[i].part->id,
parts[i].x[0], parts[i].x[1], parts[i].x[2], parts[i].v[0],
parts[i].v[1], parts[i].v[2], parts[i].a[0], parts[i].a[1],
parts[i].a[2], parts[i].mass, parts[i].dt);
found = 1;
}
if (!found) printf("## Particles[???] id=%lld not found\n", id);
}
/**
* @brief Prints the details of a given particle to stdout
*
*
* @param p The particle to print
*
*
*/
void printParticle_single ( struct part *p ) {
printf("## Particle: id=%lld, x=[%e,%e,%e], v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], h=%.3e, h_dt=%.3e, wcount=%.3e, m=%.3e, rho=%.3e, rho_dh=%.3e, div_v=%.3e, u=%.3e, dudt=%.3e, bals=%.3e, POrho2=%.3e, v_sig=%.3e, dt=%.3e\n",
p->id,
p->x[0], p->x[1], p->x[2],
p->v[0], p->v[1], p->v[2],
p->a[0], p->a[1], p->a[2],
p->h,
p->force.h_dt,
p->density.wcount,
p->mass,
p->rho, p->rho_dh,
p->density.div_v,
p->u,
p->force.u_dt,
p->force.balsara,
p->force.POrho2,
p->force.v_sig,
p->dt
);
}
void printParticle_single(struct part *p) {
printf(
"## Particle: id=%lld, x=[%e,%e,%e], v=[%.3e,%.3e,%.3e], "
"a=[%.3e,%.3e,%.3e], h=%.3e, h_dt=%.3e, wcount=%.3e, m=%.3e, rho=%.3e, "
"rho_dh=%.3e, div_v=%.3e, u=%.3e, dudt=%.3e, bals=%.3e, POrho2=%.3e, "
"v_sig=%.3e, dt=%.3e\n",
p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a[0],
p->a[1], p->a[2], p->h, p->force.h_dt, p->density.wcount, p->mass, p->rho,
p->rho_dh, p->density.div_v, p->u, p->force.u_dt, p->force.balsara,
p->force.POrho2, p->force.v_sig, p->dt);
}
src/debug.h
+ 13
− 7
View file @ 3c9a3d8a
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2012 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_DEBUG_H
#define SWIFT_DEBUG_H
/* Includes. */
#include "cell.h"
#include "part.h"
void print_cell(struct cell *c);
void printParticle(struct part *parts, long long int i, int N);
void printgParticle(struct gpart *parts, long long int i, int N);
void printParticle_single ( struct part *p );
void printParticle_single(struct part *p);
#endif /* SWIFT_DEBUG_H */
src/engine.c
+ 1955
− 1833
View file @ 3c9a3d8a
This diff is collapsed.
src/engine.h
+ 110
− 101
View file @ 3c9a3d8a
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@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_ENGINE_H
#define SWIFT_ENGINE_H
/* Some standard headers. */
#include <pthread.h>
/* Includes. */
#include "lock.h"
#include "proxy.h"
#include "runner.h"
#include "scheduler.h"
#include "space.h"
#include "task.h"
/* Some constants. */
#define engine_policy_none 0
#define engine_policy_rand 1
#define engine_policy_steal 2
#define engine_policy_keep 4
#define engine_policy_block 8
#define engine_policy_fixdt 16
#define engine_policy_multistep 32
#define engine_policy_cputight 64
#define engine_policy_mpi 128
#define engine_policy_setaffinity 256
#define engine_queue_scale 1.2
#define engine_maxtaskspercell 128
#define engine_maxproxies 64
#define engine_tasksreweight 10
#define engine_policy_none 0
#define engine_policy_rand 1
#define engine_policy_steal 2
#define engine_policy_keep 4
#define engine_policy_block 8
#define engine_policy_fixdt 16
#define engine_policy_multistep 32
#define engine_policy_cputight 64
#define engine_policy_mpi 128
#define engine_policy_setaffinity 256
#define engine_policy_paranoid 512
#define engine_queue_scale 1.2
#define engine_maxtaskspercell 128
#define engine_maxproxies 64
#define engine_tasksreweight 10
/* The rank of the engine as a global variable (for messages). */
extern int engine_rank;
/* Mini struct to link cells to density/force tasks. */
struct link {
/* The task pointer. */
struct task *t;
/* The next pointer. */
struct link *next;
};
/* The task pointer. */
struct task *t;
/* The next pointer. */
struct link *next;
};
/* Data structure for the engine. */
struct engine {
/* Number of threads on which to run. */
int nr_threads;
/* The space with which the runner is associated. */
struct space *s;
/* The runner's threads. */
struct runner *runners;
/* The running policy. */
int policy;
/* The task scheduler. */
struct scheduler sched;
/* The maximum dt to step (current). */
float dt_step;
/* The minimum dt over all particles in the system. */
float dt_min, dt_max;
/* The system time step. */
float dt, dt_orig;
/* The system energies from the previous step. */
double ekin, epot;
/* The current step number. */
int step, nullstep;
/* The number of particles updated in the previous step. */
int count_step;
/* The current system time. */
float time;
/* Data for the threads' barrier. */
pthread_mutex_t barrier_mutex;
pthread_cond_t barrier_cond;
volatile int barrier_running, barrier_launch, barrier_launchcount;
/* ID of the node this engine lives on. */
int nr_nodes, nodeID;
/* Proxies for the other nodes in this simulation. */
struct proxy *proxies;
int nr_proxies, *proxy_ind;
/* Tic at the start of a step. */
ticks tic_step;
/* Force the engine to rebuild? */
int forcerebuild, forcerepart;
/* How many steps have we done with the same set of tasks? */
int tasks_age;
/* Linked list for cell-task association. */
struct link *links;
int nr_links;
};
/* Number of threads on which to run. */
int nr_threads;
/* The space with which the runner is associated. */
struct space *s;
/* The runner's threads. */
struct runner *runners;
/* The running policy. */
int policy;
/* The task scheduler. */
struct scheduler sched;
/* The maximum dt to step (current). */
float dt_step;
/* The minimum dt over all particles in the system. */
float dt_min, dt_max;
/* The system time step. */
float dt, dt_orig;
/* The system energies from the previous step. */
double ekin, epot;
/* The current step number. */
int step, nullstep;
/* The number of particles updated in the previous step. */
int count_step;
/* The current system time. */
float time;
/* Data for the threads' barrier. */
pthread_mutex_t barrier_mutex;
pthread_cond_t barrier_cond;
volatile int barrier_running, barrier_launch, barrier_launchcount;
/* ID of the node this engine lives on. */
int nr_nodes, nodeID;
/* Proxies for the other nodes in this simulation. */
struct proxy *proxies;
int nr_proxies, *proxy_ind;
/* Tic at the start of a step. */
ticks tic_step;
/* Force the engine to rebuild? */
int forcerebuild, forcerepart;
/* How many steps have we done with the same set of tasks? */
int tasks_age;
/* Linked list for cell-task association. */
struct link *links;
int nr_links;
};
/* Function prototypes. */
void engine_barrier( struct engine *e , int tid );
void engine_init ( struct engine *e , struct space *s , float dt , int nr_threads , int nr_queues , int nr_nodes , int nodeID , int policy );
void engine_prepare ( struct engine *e );
void engine_step ( struct engine *e );
void engine_maketasks ( struct engine *e );
void engine_split ( struct engine *e , int *grid );
int engine_exchange_strays ( struct engine *e , int offset , int *ind , int N );
void engine_rebuild ( struct engine *e );
void engine_repartition ( struct engine *e );
void engine_makeproxies ( struct engine *e );
void engine_redistribute ( struct engine *e );
struct link *engine_addlink( struct engine *e , struct link *l , struct task *t );
void engine_barrier(struct engine *e, int tid);
void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
int nr_queues, int nr_nodes, int nodeID, int policy);
void engine_prepare(struct engine *e);
void engine_step(struct engine *e);
void engine_maketasks(struct engine *e);
void engine_split(struct engine *e, int *grid);
int engine_exchange_strays(struct engine *e, int offset, int *ind, int N);
void engine_rebuild(struct engine *e);
void engine_repartition(struct engine *e);
void engine_makeproxies(struct engine *e);
void engine_redistribute(struct engine *e);
struct link *engine_addlink(struct engine *e, struct link *l, struct task *t);
#endif /* SWIFT_ENGINE_H */
src/error.h
+ 30
− 11
View file @ 3c9a3d8a
......@@ -2,44 +2,63 @@
* This file is part of SWIFT.
* Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
* 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_ERROR_H
#define SWIFT_ERROR_H
/* Some standard headers. */
#include <stdio.h>
/* MPI headers. */
#ifdef WITH_MPI
#include <mpi.h>
#endif
/**
* @brief Error macro. Prints the message given in argument and aborts.
*
*/
#ifdef WITH_MPI
extern int engine_rank;
#define error(s, ...) { fprintf( stderr , "[%03i] %s:%s():%i: " s "\n" , engine_rank , __FILE__ , __FUNCTION__ , __LINE__ , ##__VA_ARGS__ ); MPI_Abort(MPI_COMM_WORLD, -1); }
extern int engine_rank;
#define error(s, ...) \
{ \
fprintf(stderr, "[%03i] %s:%s():%i: " s "\n", engine_rank, __FILE__, \
__FUNCTION__, __LINE__, ##__VA_ARGS__); \
MPI_Abort(MPI_COMM_WORLD, -1); \
}
#else
#define error(s, ...) { fprintf( stderr , "%s:%s():%i: " s "\n" , __FILE__ , __FUNCTION__ , __LINE__ , ##__VA_ARGS__ ); abort(); }
#define error(s, ...) \
{ \
fprintf(stderr, "%s:%s():%i: " s "\n", __FILE__, __FUNCTION__, __LINE__, \
##__VA_ARGS__); \
abort(); \
}
#endif
/**
* @brief Macro to print a localized message with variable arguments.
*
*/
#ifdef WITH_MPI
extern int engine_rank;
#define message(s, ...) printf( "[%03i] %s: " s "\n" , engine_rank , __FUNCTION__ , ##__VA_ARGS__ )
extern int engine_rank;
#define message(s, ...) \
printf("[%03i] %s: " s "\n", engine_rank, __FUNCTION__, ##__VA_ARGS__)
#else
#define message(s, ...) printf( "%s: " s "\n" , __FUNCTION__ , ##__VA_ARGS__ )
#define message(s, ...) printf("%s: " s "\n", __FUNCTION__, ##__VA_ARGS__)
#endif
#endif /* SWIFT_ERROR_H */
src/inline.h
+ 14
− 10
View file @ 3c9a3d8a
......@@ -2,29 +2,33 @@
* This file is part of SWIFT.
* Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
* 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_INLINE_H
#define SWIFT_INLINE_H
/**
* @brief Defines inline
* @brief Defines inline
*/
#ifndef INLINE
# if __GNUC__ && !__GNUC_STDC_INLINE__
# define INLINE extern inline
# else
# define INLINE inline
# endif
#if __GNUC__ && !__GNUC_STDC_INLINE__
#define INLINE extern inline
#else
#define INLINE inline
#endif
#endif
#endif /* SWIFT_INLINE_H */
src/kernel.h
+ 425
− 394
View file @ 3c9a3d8a
This diff is collapsed.
src/lock.h
+ 40
− 33
View file @ 3c9a3d8a
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@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_LOCK_H
#define SWIFT_LOCK_H
/* Some standard headers. */
#include <pthread.h>
/* Includes. */
#include "inline.h"
#ifdef PTHREAD_SPINLOCK
#include <pthread.h>
#define lock_type pthread_spinlock_t
#define lock_init( l ) ( pthread_spin_init( l , PTHREAD_PROCESS_PRIVATE ) != 0 )
#define lock_destroy( l ) ( pthread_spin_destroy( l ) != 0 )
#define lock_lock( l ) ( pthread_spin_lock( l ) != 0 )
#define lock_trylock( l ) ( pthread_spin_lock( l ) != 0 )
#define lock_unlock( l ) ( pthread_spin_unlock( l ) != 0 )
#define lock_unlock_blind( l ) pthread_spin_unlock( l )
#include <pthread.h>
#define lock_type pthread_spinlock_t
#define lock_init(l) (pthread_spin_init(l, PTHREAD_PROCESS_PRIVATE) != 0)
#define lock_destroy(l) (pthread_spin_destroy(l) != 0)
#define lock_lock(l) (pthread_spin_lock(l) != 0)
#define lock_trylock(l) (pthread_spin_lock(l) != 0)
#define lock_unlock(l) (pthread_spin_unlock(l) != 0)
#define lock_unlock_blind(l) pthread_spin_unlock(l)
#elif defined(PTHREAD_LOCK)
#include <pthread.h>
#define lock_type pthread_mutex_t
#define lock_init( l ) ( pthread_mutex_init( l , NULL ) != 0 )
#define lock_destroy( l ) ( pthread_mutex_destroy( l ) != 0 )
#define lock_lock( l ) ( pthread_mutex_lock( l ) != 0 )
#define lock_trylock( l ) ( pthread_mutex_trylock( l ) != 0 )
#define lock_unlock( l ) ( pthread_mutex_unlock( l ) != 0 )
#define lock_unlock_blind( l ) pthread_mutex_unlock( l )
#include <pthread.h>
#define lock_type pthread_mutex_t
#define lock_init(l) (pthread_mutex_init(l, NULL) != 0)
#define lock_destroy(l) (pthread_mutex_destroy(l) != 0)
#define lock_lock(l) (pthread_mutex_lock(l) != 0)
#define lock_trylock(l) (pthread_mutex_trylock(l) != 0)
#define lock_unlock(l) (pthread_mutex_unlock(l) != 0)
#define lock_unlock_blind(l) pthread_mutex_unlock(l)
#else
#define lock_type volatile int
#define lock_init( l ) ( *(l) = 0 )
#define lock_destroy( l ) 0
INLINE static int lock_lock ( volatile int *l ) {
while ( __sync_val_compare_and_swap( l , 0 , 1 ) != 0 );
// while( *l );
return 0;
}
#define lock_trylock( l ) ( ( *(l) ) ? 1 : __sync_val_compare_and_swap( l , 0 , 1 ) )
#define lock_unlock( l ) ( __sync_val_compare_and_swap( l , 1 , 0 ) != 1 )
#define lock_unlock_blind( l ) __sync_val_compare_and_swap( l , 1 , 0 )
#define lock_type volatile int
#define lock_init(l) (*(l) = 0)
#define lock_destroy(l) 0
INLINE static int lock_lock(volatile int *l) {
while (__sync_val_compare_and_swap(l, 0, 1) != 0)
;
// while( *l );
return 0;
}
#define lock_trylock(l) ((*(l)) ? 1 : __sync_val_compare_and_swap(l, 0, 1))
#define lock_unlock(l) (__sync_val_compare_and_swap(l, 1, 0) != 1)
#define lock_unlock_blind(l) __sync_val_compare_and_swap(l, 1, 0)
#endif
#endif /* SWIFT_LOCK_H */
src/multipole.c
+ 100
− 121
View file @ 3c9a3d8a
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@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/>.
*
*
******************************************************************************/
/* Config parameters. */
#include "../config.h"
/* Some standard headers. */
#include <float.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <limits.h>
/* MPI headers. */
#ifdef WITH_MPI
#include <mpi.h>
#include <mpi.h>
#endif
/* Local headers. */
#include "error.h"
#include "const.h"
#include "cycle.h"
#include "atomic.h"
#include "lock.h"
#include "space.h"
#include "part.h"
/* This object's header. */
#include "multipole.h"
#include "cell.h"
/**
* @brief Merge two multipoles.
......@@ -51,26 +42,23 @@
* @param ma The #multipole which will contain the merged result.
* @param mb The other #multipole.
*/
void multipole_merge ( struct multipole *ma , struct multipole *mb ) {
#if multipole_order == 1
/* Correct the position. */
float mma = ma->coeffs[0], mmb = mb->coeffs[0];
float w = 1.0f / ( mma + mmb );
for ( int k = 0 ; k < 3 ; k++ )
ma->x[k] = ( ma->x[k]*mma + mb->x[k]*mmb ) * w;
/* Add the particle to the moments. */
ma->coeffs[0] = mma + mmb;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
void multipole_merge(struct multipole *ma, struct multipole *mb) {
#if multipole_order == 1
/* Correct the position. */
float mma = ma->coeffs[0], mmb = mb->coeffs[0];
float w = 1.0f / (mma + mmb);
for (int k = 0; k < 3; k++) ma->x[k] = (ma->x[k] * mma + mb->x[k] * mmb) * w;
/* Add the particle to the moments. */
ma->coeffs[0] = mma + mmb;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
/**
* @brief Add a particle to the given multipole.
......@@ -78,26 +66,23 @@ void multipole_merge ( struct multipole *ma , struct multipole *mb ) {
* @param m The #multipole.
* @param p The #gpart.
*/
void multipole_addpart ( struct multipole *m , struct gpart *p ) {
#if multipole_order == 1
/* Correct the position. */
float mm = m->coeffs[0], mp = p->mass;
float w = 1.0f / ( mm + mp );
for ( int k = 0 ; k < 3 ; k++ )
m->x[k] = ( m->x[k]*mm + p->x[k]*mp ) * w;
/* Add the particle to the moments. */
m->coeffs[0] = mm + mp;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
void multipole_addpart(struct multipole *m, struct gpart *p) {
#if multipole_order == 1
/* Correct the position. */
float mm = m->coeffs[0], mp = p->mass;
float w = 1.0f / (mm + mp);
for (int k = 0; k < 3; k++) m->x[k] = (m->x[k] * mm + p->x[k] * mp) * w;
/* Add the particle to the moments. */
m->coeffs[0] = mm + mp;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
/**
* @brief Add a group of particles to the given multipole.
......@@ -106,37 +91,34 @@ void multipole_addpart ( struct multipole *m , struct gpart *p ) {
* @param p The #gpart array.
* @param N Number of parts to add.
*/
void multipole_addparts ( struct multipole *m , struct gpart *p , int N ) {
#if multipole_order == 1
/* Get the combined mass and positions. */
double xp[3] = { 0.0 , 0.0 , 0.0 };
float mp = 0.0f, w;
for ( int k = 0 ; k < N ; k++ ) {
w = p[k].mass;
mp += w;
xp[0] += p[k].x[0] * w;
xp[1] += p[k].x[1] * w;
xp[2] += p[k].x[2] * w;
}
/* Correct the position. */
float mm = m->coeffs[0];
w = 1.0f / ( mm + mp );
for ( int k = 0 ; k < 3 ; k++ )
m->x[k] = ( m->x[k]*mm + xp[k] ) * w;
/* Add the particle to the moments. */
m->coeffs[0] = mm + mp;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
void multipole_addparts(struct multipole *m, struct gpart *p, int N) {
#if multipole_order == 1
/* Get the combined mass and positions. */
double xp[3] = {0.0, 0.0, 0.0};
float mp = 0.0f, w;
for (int k = 0; k < N; k++) {
w = p[k].mass;
mp += w;
xp[0] += p[k].x[0] * w;
xp[1] += p[k].x[1] * w;
xp[2] += p[k].x[2] * w;
}
/* Correct the position. */
float mm = m->coeffs[0];
w = 1.0f / (mm + mp);
for (int k = 0; k < 3; k++) m->x[k] = (m->x[k] * mm + xp[k]) * w;
/* Add the particle to the moments. */
m->coeffs[0] = mm + mp;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
/**
* @brief Init a multipole from a set of particles.
......@@ -145,46 +127,43 @@ void multipole_addparts ( struct multipole *m , struct gpart *p , int N ) {
* @param parts The #gpart.
* @param N The number of particles.
*/
void multipole_init ( struct multipole *m , struct gpart *parts , int N ) {
#if multipole_order == 1
float mass = 0.0f, w;
double x[3] = { 0.0 , 0.0 , 0.0 };
int k;
/* Collect the particle data. */
for ( k = 0 ; k < N ; k++ ) {
w = parts[k].mass;
mass += w;
x[0] += parts[k].x[0] * w;
x[1] += parts[k].x[1] * w;
x[2] += parts[k].x[2] * w;
}
/* Store the data on the multipole. */
m->coeffs[0] = mass;
m->x[0] = x[0] / mass;
m->x[1] = x[1] / mass;
m->x[2] = x[2] / mass;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
void multipole_init(struct multipole *m, struct gpart *parts, int N) {
#if multipole_order == 1
float mass = 0.0f, w;
double x[3] = {0.0, 0.0, 0.0};
int k;
/* Collect the particle data. */
for (k = 0; k < N; k++) {
w = parts[k].mass;
mass += w;
x[0] += parts[k].x[0] * w;
x[1] += parts[k].x[1] * w;
x[2] += parts[k].x[2] * w;
}
/* Store the data on the multipole. */
m->coeffs[0] = mass;
m->x[0] = x[0] / mass;
m->x[1] = x[1] / mass;
m->x[2] = x[2] / mass;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
/**
* @brief Reset the data of a #multipole.
*
* @param m The #multipole.
*/
void multipole_reset ( struct multipole *m ) {
/* Just bzero the struct. */
bzero( m , sizeof(struct multipole) );
}
void multipole_reset(struct multipole *m) {
/* Just bzero the struct. */
bzero(m, sizeof(struct multipole));
}
src/multipole.h
+ 96
− 93
View file @ 3c9a3d8a
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@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_MULTIPOLE_H
#define SWIFT_MULTIPOLE_H
/* Some constants. */
#define multipole_order 1
/* Some standard headers. */
#include <math.h>
/* Includes. */
#include "const.h"
#include "inline.h"
#include "kernel.h"
#include "part.h"
/* Some constants. */
#define multipole_order 1
/* Multipole struct. */
struct multipole {
/* Multipole location. */
double x[3];
/* Acceleration on this multipole. */
float a[3];
/* Multipole coefficients. */
float coeffs[ multipole_order*multipole_order ];
};
/* Multipole function prototypes. */
static void multipole_iact_mm ( struct multipole *ma , struct multipole *mb , double *shift );
void multipole_merge ( struct multipole *ma , struct multipole *mb );
void multipole_addpart ( struct multipole *m , struct gpart *p );
void multipole_addparts ( struct multipole *m , struct gpart *p , int N );
void multipole_init ( struct multipole *m , struct gpart *parts , int N );
void multipole_reset ( struct multipole *m );
/* Multipole location. */
double x[3];
/* Acceleration on this multipole. */
float a[3];
/* Multipole coefficients. */
float coeffs[multipole_order * multipole_order];
};
/* Multipole function prototypes. */
static void multipole_iact_mm(struct multipole *ma, struct multipole *mb,
double *shift);
void multipole_merge(struct multipole *ma, struct multipole *mb);
void multipole_addpart(struct multipole *m, struct gpart *p);
void multipole_addparts(struct multipole *m, struct gpart *p, int N);
void multipole_init(struct multipole *m, struct gpart *parts, int N);
void multipole_reset(struct multipole *m);
#include <math.h>
#include "kernel.h"
/**
* @brief Compute the pairwise interaction between two multipoles.
*
......@@ -55,41 +60,40 @@ void multipole_reset ( struct multipole *m );
* @param mb The second #multipole.
* @param shift The periodicity correction.
*/
__attribute__ ((always_inline)) INLINE static void multipole_iact_mm ( struct multipole *ma , struct multipole *mb , double *shift ) {
float dx[3], ir, r, r2 = 0.0f, acc;
int k;
/* Compute the multipole distance. */
for ( k = 0 ; k < 3 ; k++ ) {
dx[k] = ma->x[k] - mb->x[k] - shift[k];
r2 += dx[k]*dx[k];
}
/* Compute the normalized distance vector. */
ir = 1.0f / sqrtf( r2 );
r = r2 * ir;
/* Evaluate the gravity kernel. */
kernel_grav_eval( r , &acc );
/* Scale the acceleration. */
acc *= const_G * ir * ir * ir;
/* Compute the forces on both multipoles. */
#if multipole_order == 1
float mma = ma->coeffs[0], mmb = mb->coeffs[0];
for ( k = 0 ; k < 3 ; k++ ) {
ma->a[k] -= dx[k] * acc * mmb;
mb->a[k] += dx[k] * acc * mma;
}
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
__attribute__((always_inline)) INLINE static void multipole_iact_mm(
struct multipole *ma, struct multipole *mb, double *shift) {
float dx[3], ir, r, r2 = 0.0f, acc;
int k;
/* Compute the multipole distance. */
for (k = 0; k < 3; k++) {
dx[k] = ma->x[k] - mb->x[k] - shift[k];
r2 += dx[k] * dx[k];
}
/* Compute the normalized distance vector. */
ir = 1.0f / sqrtf(r2);
r = r2 * ir;
/* Evaluate the gravity kernel. */
kernel_grav_eval(r, &acc);
/* Scale the acceleration. */
acc *= const_G * ir * ir * ir;
/* Compute the forces on both multipoles. */
#if multipole_order == 1
float mma = ma->coeffs[0], mmb = mb->coeffs[0];
for (k = 0; k < 3; k++) {
ma->a[k] -= dx[k] * acc * mmb;
mb->a[k] += dx[k] * acc * mma;
}
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
/**
* @brief Compute the interaction of a multipole on a particle.
......@@ -98,36 +102,35 @@ __attribute__ ((always_inline)) INLINE static void multipole_iact_mm ( struct mu
* @param p The #gpart.
* @param shift The periodicity correction.
*/
__attribute__ ((always_inline)) INLINE static void multipole_iact_mp ( struct multipole *m , struct gpart *p , double *shift ) {
float dx[3], ir, r, r2 = 0.0f, acc;
int k;
/* Compute the multipole distance. */
for ( k = 0 ; k < 3 ; k++ ) {
dx[k] = m->x[k] - p->x[k] - shift[k];
r2 += dx[k]*dx[k];
}
/* Compute the normalized distance vector. */
ir = 1.0f / sqrtf( r2 );
r = r2 * ir;
/* Evaluate the gravity kernel. */
kernel_grav_eval( r , &acc );
/* Scale the acceleration. */
acc *= const_G * ir * ir * ir * m->coeffs[0];
/* Compute the forces on both multipoles. */
#if multipole_order == 1
for ( k = 0 ; k < 3 ; k++ )
p->a[k] += dx[k] * acc;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
__attribute__((always_inline)) INLINE static void multipole_iact_mp(
struct multipole *m, struct gpart *p, double *shift) {
float dx[3], ir, r, r2 = 0.0f, acc;
int k;
/* Compute the multipole distance. */
for (k = 0; k < 3; k++) {
dx[k] = m->x[k] - p->x[k] - shift[k];
r2 += dx[k] * dx[k];
}
/* Compute the normalized distance vector. */
ir = 1.0f / sqrtf(r2);
r = r2 * ir;
/* Evaluate the gravity kernel. */
kernel_grav_eval(r, &acc);
/* Scale the acceleration. */
acc *= const_G * ir * ir * ir * m->coeffs[0];
/* Compute the forces on both multipoles. */
#if multipole_order == 1
for (k = 0; k < 3; k++) p->a[k] += dx[k] * acc;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
#endif /* SWIFT_MULTIPOLE_H */
src/parallel_io.c
+ 238
− 228
View file @ 3c9a3d8a
This diff is collapsed.
src/parallel_io.h
+ 22
− 6
View file @ 3c9a3d8a
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment