Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
SWIFT
SWIFTsim
Commits
0167dd62
Commit
0167dd62
authored
May 10, 2017
by
James Willis
Browse files
Merge branch 'master' into doself2-vectorisation
Conflicts: src/cache.h src/runner_doiact_vec.c
parents
dacc84d8
0602056b
Changes
64
Expand all
Hide whitespace changes
Inline
Side-by-side
.gitignore
View file @
0167dd62
...
...
@@ -91,6 +91,7 @@ theory/SPH/*.pdf
theory/paper_pasc/pasc_paper.pdf
theory/Multipoles/fmm.pdf
theory/Multipoles/fmm_standalone.pdf
theory/Multipoles/potential.pdf
m4/libtool.m4
m4/ltoptions.m4
...
...
README
View file @
0167dd62
...
...
@@ -27,6 +27,7 @@ Valid options are:
-f {int} Overwrite the CPU frequency (Hz) to be used for time measurements.
-g Run with an external gravitational potential.
-G Run with self-gravity.
-M Reconstruct the multipoles every time-step.
-n {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop.
-s Run with hydrodynamics.
-S Run with stars.
...
...
configure.ac
View file @
0167dd62
...
...
@@ -189,6 +189,18 @@ if test "$enable_task_debugging" = "yes"; then
AC_DEFINE([SWIFT_DEBUG_TASKS],1,[Enable task debugging])
fi
# Check if the general timers are switched on.
AC_ARG_ENABLE([timers],
[AS_HELP_STRING([--enable-timers],
[Activate the basic timers @<:@yes/no@:>@]
)],
[enable_timers="$enableval"],
[enable_timers="no"]
)
if test "$enable_timers" = "yes"; then
AC_DEFINE([SWIFT_USE_TIMERS],1,[Enable individual timers])
fi
# Check if expensive debugging is on.
AC_ARG_ENABLE([debugging-checks],
[AS_HELP_STRING([--enable-debugging-checks],
...
...
@@ -215,6 +227,21 @@ elif test "$gravity_force_checks" != "no"; then
AC_DEFINE_UNQUOTED([SWIFT_GRAVITY_FORCE_CHECKS], [$enableval] ,[Enable gravity brute-force checks])
fi
# Check if we want to zero the gravity forces for all particles below some ID.
AC_ARG_ENABLE([no-gravity-below-id],
[AS_HELP_STRING([--enable-no-gravity-below-id],
[Zeros the gravitational acceleration of all particles with an ID smaller than @<:@N@:>@]
)],
[no_gravity_below_id="$enableval"],
[no_gravity_below_id="no"]
)
if test "$no_gravity_below_id" == "yes"; then
AC_MSG_ERROR(Need to specify the ID below which particles get zero forces when using --enable-no-gravity-below-id!)
elif test "$no_gravity_below_id" != "no"; then
AC_DEFINE_UNQUOTED([SWIFT_NO_GRAVITY_BELOW_ID], [$enableval] ,[Particles with smaller ID than this will have zero gravity forces])
fi
# Define HAVE_POSIX_MEMALIGN if it works.
AX_FUNC_POSIX_MEMALIGN
...
...
@@ -854,12 +881,15 @@ AC_MSG_RESULT([
Adiabatic index : $with_gamma
Riemann solver : $with_riemann
Cooling function : $with_cooling
External potential : $with_potential
Multipole order : $with_multipole_order
Task debugging : $enable_task_debugging
Debugging checks : $enable_debugging_checks
Gravity checks : $gravity_force_checks
External potential : $with_potential
Multipole order : $with_multipole_order
No gravity below ID : $no_gravity_below_id
Individual timers : $enable_timers
Task debugging : $enable_task_debugging
Debugging checks : $enable_debugging_checks
Gravity checks : $gravity_force_checks
])
# Make sure the latest git revision string gets included
...
...
examples/Makefile.am
View file @
0167dd62
...
...
@@ -16,7 +16,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Common flags
MYFLAGS
=
-DTIMER
MYFLAGS
=
# Add the source directory and debug to CFLAGS
AM_CFLAGS
=
-I
$(top_srcdir)
/src
$(HDF5_CPPFLAGS)
...
...
examples/main.c
View file @
0167dd62
...
...
@@ -82,6 +82,8 @@ void print_help_message() {
"Run with an external gravitational potential."
);
printf
(
" %2s %8s %s
\n
"
,
"-F"
,
""
,
"Run with feedback."
);
printf
(
" %2s %8s %s
\n
"
,
"-G"
,
""
,
"Run with self-gravity."
);
printf
(
" %2s %8s %s
\n
"
,
"-M"
,
""
,
"Reconstruct the multipoles every time-step."
);
printf
(
" %2s %8s %s
\n
"
,
"-n"
,
"{int}"
,
"Execute a fixed number of time steps. When unset use the time_end "
"parameter to stop."
);
...
...
@@ -164,6 +166,7 @@ int main(int argc, char *argv[]) {
int
with_stars
=
0
;
int
with_fp_exceptions
=
0
;
int
with_drift_all
=
0
;
int
with_mpole_reconstruction
=
0
;
int
verbose
=
0
;
int
nr_threads
=
1
;
int
with_verbose_timers
=
0
;
...
...
@@ -172,7 +175,8 @@ int main(int argc, char *argv[]) {
/* Parse the parameters */
int
c
;
while
((
c
=
getopt
(
argc
,
argv
,
"acCdDef:FgGhn:sSt:Tv:y:"
))
!=
-
1
)
switch
(
c
)
{
while
((
c
=
getopt
(
argc
,
argv
,
"acCdDef:FgGhMn:sSt:Tv:y:"
))
!=
-
1
)
switch
(
c
)
{
case
'a'
:
with_aff
=
1
;
break
;
...
...
@@ -210,6 +214,9 @@ int main(int argc, char *argv[]) {
case
'h'
:
if
(
myrank
==
0
)
print_help_message
();
return
0
;
case
'M'
:
with_mpole_reconstruction
=
1
;
break
;
case
'n'
:
if
(
sscanf
(
optarg
,
"%d"
,
&
nsteps
)
!=
1
)
{
if
(
myrank
==
0
)
printf
(
"Error parsing fixed number of steps.
\n
"
);
...
...
@@ -521,6 +528,8 @@ int main(int argc, char *argv[]) {
/* Construct the engine policy */
int
engine_policies
=
ENGINE_POLICY
|
engine_policy_steal
;
if
(
with_drift_all
)
engine_policies
|=
engine_policy_drift_all
;
if
(
with_mpole_reconstruction
)
engine_policies
|=
engine_policy_reconstruct_mpoles
;
if
(
with_hydro
)
engine_policies
|=
engine_policy_hydro
;
if
(
with_self_gravity
)
engine_policies
|=
engine_policy_self_gravity
;
if
(
with_external_gravity
)
engine_policies
|=
engine_policy_external_gravity
;
...
...
@@ -601,35 +610,25 @@ int main(int argc, char *argv[]) {
engine_dump_snapshot
(
&
e
);
/* Legend */
if
(
myrank
==
0
)
{
if
(
myrank
==
0
)
printf
(
"# %6s %14s %14s %10s %10s %10s %16s [%s]
\n
"
,
"Step"
,
"Time"
,
"Time-step"
,
"Updates"
,
"g-Updates"
,
"s-Updates"
,
"Wall-clock time"
,
clocks_getunit
());
if
(
with_verbose_timers
)
{
printf
(
"timers: "
);
for
(
int
k
=
0
;
k
<
timer_count
;
k
++
)
printf
(
"%s
\t
"
,
timers_names
[
k
]);
printf
(
"
\n
"
);
}
}
/* File for the timers */
if
(
with_verbose_timers
)
timers_open_file
(
myrank
);
/* Main simulation loop */
for
(
int
j
=
0
;
!
engine_is_done
(
&
e
)
&&
e
.
step
-
1
!=
nsteps
;
j
++
)
{
/* Reset timers */
timers_reset
(
timers_mask
_all
);
timers_reset_all
(
);
/* Take a step. */
engine_step
(
&
e
);
/* Print the timers. */
if
(
with_verbose_timers
)
{
printf
(
"timers: "
);
for
(
int
k
=
0
;
k
<
timer_count
;
k
++
)
printf
(
"%.3f
\t
"
,
clocks_from_ticks
(
timers
[
k
]));
printf
(
"
\n
"
);
timers_reset
(
0xFFFFFFFFllu
);
}
if
(
with_verbose_timers
)
timers_print
(
e
.
step
);
#ifdef SWIFT_DEBUG_TASKS
/* Dump the task data using the given frequency. */
...
...
@@ -735,6 +734,7 @@ int main(int argc, char *argv[]) {
#endif
/* Clean everything */
if
(
with_verbose_timers
)
timers_close_file
();
engine_clean
(
&
e
);
free
(
params
);
...
...
src/Makefile.am
View file @
0167dd62
...
...
@@ -16,7 +16,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Add the debug flag to the whole thing
AM_CFLAGS
=
-DTIMER
$(HDF5_CPPFLAGS)
AM_CFLAGS
=
$(HDF5_CPPFLAGS)
# Assign a "safe" version number
AM_LDFLAGS
=
$(HDF5_LDFLAGS)
$(FFTW_LIBS)
-version-info
0:0:0
...
...
@@ -46,7 +46,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h
\
sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h
\
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h
\
vector_power.h collectgroup.h hydro_space.h sort_part.h
gravity_softened_derivatives.h
vector_power.h collectgroup.h hydro_space.h sort_part.h
# Common source files
AM_SOURCES
=
space.c runner.c queue.c task.c cell.c engine.c
\
...
...
src/atomic.h
View file @
0167dd62
...
...
@@ -19,6 +19,9 @@
#ifndef SWIFT_ATOMIC_H
#define SWIFT_ATOMIC_H
/* Config parameters. */
#include
"../config.h"
/* Includes. */
#include
"inline.h"
...
...
src/cache.h
View file @
0167dd62
...
...
@@ -34,6 +34,7 @@
#define C2_CACHE_SIZE (NUM_VEC_PROC * VEC_SIZE * 6) + (NUM_VEC_PROC * VEC_SIZE)
#define C2_CACHE_ALIGN sizeof(float) * VEC_SIZE
#ifdef WITH_VECTORIZATION
/* Cache struct to hold a local copy of a cells' particle
* properties required for density/force calculations.*/
struct
cache
{
...
...
@@ -178,10 +179,10 @@ __attribute__((always_inline)) INLINE void cache_read_particles(
#if defined(GADGET2_SPH)
/* Shift the particles positions to a local frame so single precision can be
* used instead of double precision. */
/* Shift the particles positions to a local frame so single precision can be
* used instead of double precision. */
#if defined(WITH_VECTORIZATION) && defined(__ICC)
#pragma
sim
d
#pragma
vector aligne
d
#endif
for
(
int
i
=
0
;
i
<
ci
->
count
;
i
++
)
{
ci_cache
->
x
[
i
]
=
ci
->
parts
[
i
].
x
[
0
]
-
ci
->
loc
[
0
];
...
...
@@ -380,7 +381,7 @@ __attribute__((always_inline)) INLINE void cache_read_two_partial_cells_sorted(
* used instead of double precision. Also shift the cell ci, particles positions
* due to BCs but leave cell cj. */
#if defined(WITH_VECTORIZATION) && defined(__ICC)
#pragma
sim
d
#pragma
vector aligne
d
#endif
for
(
int
i
=
first_pi_align
;
i
<
ci
->
count
;
i
++
)
{
/* Make sure ci_cache is filled from the first element. */
...
...
@@ -397,14 +398,24 @@ __attribute__((always_inline)) INLINE void cache_read_two_partial_cells_sorted(
ci_cache
->
vz
[
ci_cache_idx
]
=
ci
->
parts
[
idx
].
v
[
2
];
}
/* Pad cache with fake particles that exist outside the cell so will not interact.*/
float
fake_pix
=
2
.
0
f
*
ci_cache
->
x
[
ci
->
count
-
1
];
/* Pad cache with fake particles that exist outside the cell so will not
* interact.*/
float
fake_pix
=
2
.
0
f
*
ci
->
parts
[
sort_i
[
ci
->
count
-
1
].
i
].
x
[
0
];
for
(
int
i
=
ci
->
count
-
first_pi_align
;
i
<
ci
->
count
-
first_pi_align
+
VEC_SIZE
;
i
++
)
i
<
ci
->
count
-
first_pi_align
+
VEC_SIZE
;
i
++
)
{
ci_cache
->
x
[
i
]
=
fake_pix
;
ci_cache
->
y
[
i
]
=
1
.
f
;
ci_cache
->
z
[
i
]
=
1
.
f
;
ci_cache
->
h
[
i
]
=
1
.
f
;
ci_cache
->
m
[
i
]
=
1
.
f
;
ci_cache
->
vx
[
i
]
=
1
.
f
;
ci_cache
->
vy
[
i
]
=
1
.
f
;
ci_cache
->
vz
[
i
]
=
1
.
f
;
}
#if defined(WITH_VECTORIZATION) && defined(__ICC)
#pragma
sim
d
#pragma
vector aligne
d
#endif
for
(
int
i
=
0
;
i
<=
last_pj_align
;
i
++
)
{
idx
=
sort_j
[
i
].
i
;
...
...
@@ -419,10 +430,20 @@ __attribute__((always_inline)) INLINE void cache_read_two_partial_cells_sorted(
cj_cache
->
vz
[
i
]
=
cj
->
parts
[
idx
].
v
[
2
];
}
/* Pad cache with fake particles that exist outside the cell so will not interact.*/
float
fake_pjx
=
2
.
0
f
*
cj_cache
->
x
[
last_pj_align
];
for
(
int
i
=
last_pj_align
+
1
;
i
<
last_pj_align
+
1
+
VEC_SIZE
;
i
++
)
/* Pad cache with fake particles that exist outside the cell so will not
* interact.*/
float
fake_pjx
=
2
.
0
f
*
cj
->
parts
[
sort_j
[
cj
->
count
-
1
].
i
].
x
[
0
];
for
(
int
i
=
last_pj_align
+
1
;
i
<
last_pj_align
+
1
+
VEC_SIZE
;
i
++
)
{
cj_cache
->
x
[
i
]
=
fake_pjx
;
cj_cache
->
y
[
i
]
=
1
.
f
;
cj_cache
->
z
[
i
]
=
1
.
f
;
cj_cache
->
h
[
i
]
=
1
.
f
;
cj_cache
->
m
[
i
]
=
1
.
f
;
cj_cache
->
vx
[
i
]
=
1
.
f
;
cj_cache
->
vy
[
i
]
=
1
.
f
;
cj_cache
->
vz
[
i
]
=
1
.
f
;
}
}
/* @brief Clean the memory allocated by a #cache object.
...
...
@@ -443,4 +464,6 @@ static INLINE void cache_clean(struct cache *c) {
}
}
#endif
/* WITH_VECTORIZATION */
#endif
/* SWIFT_CACHE_H */
src/cell.c
View file @
0167dd62
...
...
@@ -129,6 +129,7 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
temp
->
depth
=
c
->
depth
+
1
;
temp
->
split
=
0
;
temp
->
dx_max
=
0
.
f
;
temp
->
dx_max_sort
=
0
.
f
;
temp
->
nodeID
=
c
->
nodeID
;
temp
->
parent
=
c
;
c
->
progeny
[
k
]
=
temp
;
...
...
@@ -1103,33 +1104,93 @@ void cell_reset_task_counters(struct cell *c) {
}
/**
* @brief Checks whether the cells are direct neighbours ot not. Both cells have
* to be of the same size
* @brief Recursively construct all the multipoles in a cell hierarchy.
*
* @param ci First #cell.
* @param cj Second #cell.
*
* @todo Deal with periodicity.
* @param c The #cell.
* @param ti_current The current integer time.
*/
int
cell_are_neighbours
(
const
struct
cell
*
restrict
ci
,
const
struct
cell
*
restrict
cj
)
{
void
cell_make_multipoles
(
struct
cell
*
c
,
integertime_t
ti_current
)
{
#ifdef SWIFT_DEBUG_CHECKS
if
(
ci
->
width
[
0
]
!=
cj
->
width
[
0
])
error
(
"Cells of different size !"
);
#endif
/* Reset everything */
gravity_reset
(
c
->
multipole
);
/* Maximum allowed distance */
const
double
min_dist
=
1
.
2
*
ci
->
width
[
0
];
/* 1.2 accounts for rounding errors */
if
(
c
->
split
)
{
/* (Manhattan) Distance between the cells */
for
(
int
k
=
0
;
k
<
3
;
k
++
)
{
const
double
center_i
=
ci
->
loc
[
k
];
const
double
center_j
=
cj
->
loc
[
k
];
if
(
fabs
(
center_i
-
center_j
)
>
min_dist
)
return
0
;
/* Compute CoM of all progenies */
double
CoM
[
3
]
=
{
0
.,
0
.,
0
.};
double
mass
=
0
.;
for
(
int
k
=
0
;
k
<
8
;
++
k
)
{
if
(
c
->
progeny
[
k
]
!=
NULL
)
{
const
struct
gravity_tensors
*
m
=
c
->
progeny
[
k
]
->
multipole
;
CoM
[
0
]
+=
m
->
CoM
[
0
]
*
m
->
m_pole
.
M_000
;
CoM
[
1
]
+=
m
->
CoM
[
1
]
*
m
->
m_pole
.
M_000
;
CoM
[
2
]
+=
m
->
CoM
[
2
]
*
m
->
m_pole
.
M_000
;
mass
+=
m
->
m_pole
.
M_000
;
}
}
c
->
multipole
->
CoM
[
0
]
=
CoM
[
0
]
/
mass
;
c
->
multipole
->
CoM
[
1
]
=
CoM
[
1
]
/
mass
;
c
->
multipole
->
CoM
[
2
]
=
CoM
[
2
]
/
mass
;
/* Now shift progeny multipoles and add them up */
struct
multipole
temp
;
double
r_max
=
0
.;
for
(
int
k
=
0
;
k
<
8
;
++
k
)
{
if
(
c
->
progeny
[
k
]
!=
NULL
)
{
const
struct
cell
*
cp
=
c
->
progeny
[
k
];
const
struct
multipole
*
m
=
&
cp
->
multipole
->
m_pole
;
/* Contribution to multipole */
gravity_M2M
(
&
temp
,
m
,
c
->
multipole
->
CoM
,
cp
->
multipole
->
CoM
);
gravity_multipole_add
(
&
c
->
multipole
->
m_pole
,
&
temp
);
/* Upper limit of max CoM<->gpart distance */
const
double
dx
=
c
->
multipole
->
CoM
[
0
]
-
cp
->
multipole
->
CoM
[
0
];
const
double
dy
=
c
->
multipole
->
CoM
[
1
]
-
cp
->
multipole
->
CoM
[
1
];
const
double
dz
=
c
->
multipole
->
CoM
[
2
]
-
cp
->
multipole
->
CoM
[
2
];
const
double
r2
=
dx
*
dx
+
dy
*
dy
+
dz
*
dz
;
r_max
=
max
(
r_max
,
cp
->
multipole
->
r_max
+
sqrt
(
r2
));
}
}
/* Alternative upper limit of max CoM<->gpart distance */
const
double
dx
=
c
->
multipole
->
CoM
[
0
]
>
c
->
loc
[
0
]
+
c
->
width
[
0
]
/
2
.
?
c
->
multipole
->
CoM
[
0
]
-
c
->
loc
[
0
]
:
c
->
loc
[
0
]
+
c
->
width
[
0
]
-
c
->
multipole
->
CoM
[
0
];
const
double
dy
=
c
->
multipole
->
CoM
[
1
]
>
c
->
loc
[
1
]
+
c
->
width
[
1
]
/
2
.
?
c
->
multipole
->
CoM
[
1
]
-
c
->
loc
[
1
]
:
c
->
loc
[
1
]
+
c
->
width
[
1
]
-
c
->
multipole
->
CoM
[
1
];
const
double
dz
=
c
->
multipole
->
CoM
[
2
]
>
c
->
loc
[
2
]
+
c
->
width
[
2
]
/
2
.
?
c
->
multipole
->
CoM
[
2
]
-
c
->
loc
[
2
]
:
c
->
loc
[
2
]
+
c
->
width
[
2
]
-
c
->
multipole
->
CoM
[
2
];
/* Take minimum of both limits */
c
->
multipole
->
r_max
=
min
(
r_max
,
sqrt
(
dx
*
dx
+
dy
*
dy
+
dz
*
dz
));
}
else
{
if
(
c
->
gcount
>
0
)
{
gravity_P2M
(
c
->
multipole
,
c
->
gparts
,
c
->
gcount
);
const
double
dx
=
c
->
multipole
->
CoM
[
0
]
>
c
->
loc
[
0
]
+
c
->
width
[
0
]
/
2
.
?
c
->
multipole
->
CoM
[
0
]
-
c
->
loc
[
0
]
:
c
->
loc
[
0
]
+
c
->
width
[
0
]
-
c
->
multipole
->
CoM
[
0
];
const
double
dy
=
c
->
multipole
->
CoM
[
1
]
>
c
->
loc
[
1
]
+
c
->
width
[
1
]
/
2
.
?
c
->
multipole
->
CoM
[
1
]
-
c
->
loc
[
1
]
:
c
->
loc
[
1
]
+
c
->
width
[
1
]
-
c
->
multipole
->
CoM
[
1
];
const
double
dz
=
c
->
multipole
->
CoM
[
2
]
>
c
->
loc
[
2
]
+
c
->
width
[
2
]
/
2
.
?
c
->
multipole
->
CoM
[
2
]
-
c
->
loc
[
2
]
:
c
->
loc
[
2
]
+
c
->
width
[
2
]
-
c
->
multipole
->
CoM
[
2
];
c
->
multipole
->
r_max
=
sqrt
(
dx
*
dx
+
dy
*
dy
+
dz
*
dz
);
}
else
{
gravity_multipole_init
(
&
c
->
multipole
->
m_pole
);
c
->
multipole
->
CoM
[
0
]
=
c
->
loc
[
0
]
+
c
->
width
[
0
]
/
2
.;
c
->
multipole
->
CoM
[
1
]
=
c
->
loc
[
1
]
+
c
->
width
[
1
]
/
2
.;
c
->
multipole
->
CoM
[
2
]
=
c
->
loc
[
2
]
+
c
->
width
[
2
]
/
2
.;
c
->
multipole
->
r_max
=
0
.;
}
}
return
1
;
c
->
ti_old_multipole
=
ti_current
;
}
/**
...
...
@@ -1145,6 +1206,8 @@ void cell_check_multipole(struct cell *c, void *data) {
struct
gravity_tensors
ma
;
const
double
tolerance
=
1e-3
;
/* Relative */
return
;
/* First recurse */
if
(
c
->
split
)
for
(
int
k
=
0
;
k
<
8
;
k
++
)
...
...
@@ -1244,28 +1307,44 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
/* Un-skip the density tasks involved with this cell. */
for
(
struct
link
*
l
=
c
->
density
;
l
!=
NULL
;
l
=
l
->
next
)
{
struct
task
*
t
=
l
->
t
;
const
struct
cell
*
ci
=
t
->
ci
;
const
struct
cell
*
cj
=
t
->
cj
;
struct
cell
*
ci
=
t
->
ci
;
struct
cell
*
cj
=
t
->
cj
;
scheduler_activate
(
s
,
t
);
/* Set the correct sorting flags */
if
(
t
->
type
==
task_type_pair
)
{
if
(
ci
->
dx_max_sort
>
space_maxreldx
*
ci
->
dmin
)
{
for
(
struct
cell
*
finger
=
ci
;
finger
!=
NULL
;
finger
=
finger
->
parent
)
finger
->
sorted
=
0
;
}
if
(
cj
->
dx_max_sort
>
space_maxreldx
*
cj
->
dmin
)
{
for
(
struct
cell
*
finger
=
cj
;
finger
!=
NULL
;
finger
=
finger
->
parent
)
finger
->
sorted
=
0
;
}
if
(
!
(
ci
->
sorted
&
(
1
<<
t
->
flags
)))
{
atomic_or
(
&
ci
->
sorts
->
flags
,
(
1
<<
t
->
flags
));
#ifdef SWIFT_DEBUG_CHECKS
if
(
!
(
ci
->
sorts
->
flags
&
(
1
<<
t
->
flags
)))
error
(
"bad flags in sort task."
);
#endif
scheduler_activate
(
s
,
ci
->
sorts
);
if
(
ci
->
nodeID
==
engine_rank
)
scheduler_activate
(
s
,
ci
->
drift
);
}
if
(
!
(
cj
->
sorted
&
(
1
<<
t
->
flags
)))
{
atomic_or
(
&
cj
->
sorts
->
flags
,
(
1
<<
t
->
flags
));
#ifdef SWIFT_DEBUG_CHECKS
if
(
!
(
cj
->
sorts
->
flags
&
(
1
<<
t
->
flags
)))
error
(
"bad flags in sort task."
);
#endif
scheduler_activate
(
s
,
cj
->
sorts
);
if
(
cj
->
nodeID
==
engine_rank
)
scheduler_activate
(
s
,
cj
->
drift
);
}
}
/*
Check whether there was too much particle motion
*/
/*
Only interested in pair interactions as of here.
*/
if
(
t
->
type
==
task_type_pair
||
t
->
type
==
task_type_sub_pair
)
{
if
(
t
->
tight
&&
(
max
(
ci
->
h_max
,
cj
->
h_max
)
+
ci
->
dx_max
+
cj
->
dx_max
>
cj
->
dmin
||
ci
->
dx_max
>
space_maxreldx
*
ci
->
h_max
||
c
j
->
dx_max
>
space_maxreldx
*
cj
->
h_max
)
)
/* Check whether there was too much particle motion, i.e. the
cell neighbour conditions were violated. */
if
(
max
(
ci
->
h_max
,
cj
->
h_max
)
+
c
i
->
dx_max
+
cj
->
dx_max
>
cj
->
dmin
)
rebuild
=
1
;
#ifdef WITH_MPI
...
...
@@ -1287,10 +1366,12 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if
(
l
==
NULL
)
error
(
"Missing link to send_xv task."
);
scheduler_activate
(
s
,
l
->
t
);
if
(
cj
->
super
->
drift
)
scheduler_activate
(
s
,
cj
->
super
->
drift
);
/* Drift both cells, the foreign one at the level which it is sent. */
if
(
l
->
t
->
ci
->
drift
)
scheduler_activate
(
s
,
l
->
t
->
ci
->
drift
);
else
error
(
"Drift task missing !"
);
if
(
t
->
type
==
task_type_pair
)
scheduler_activate
(
s
,
cj
->
drift
);
if
(
cell_is_active
(
cj
,
e
))
{
for
(
l
=
cj
->
send_rho
;
l
!=
NULL
&&
l
->
t
->
cj
->
nodeID
!=
ci
->
nodeID
;
...
...
@@ -1323,10 +1404,12 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if
(
l
==
NULL
)
error
(
"Missing link to send_xv task."
);
scheduler_activate
(
s
,
l
->
t
);
if
(
ci
->
super
->
drift
)
scheduler_activate
(
s
,
ci
->
super
->
drift
);
/* Drift both cells, the foreign one at the level which it is sent. */
if
(
l
->
t
->
ci
->
drift
)
scheduler_activate
(
s
,
l
->
t
->
ci
->
drift
);
else
error
(
"Drift task missing !"
);
if
(
t
->
type
==
task_type_pair
)
scheduler_activate
(
s
,
ci
->
drift
);
if
(
cell_is_active
(
ci
,
e
))
{
for
(
l
=
ci
->
send_rho
;
l
!=
NULL
&&
l
->
t
->
cj
->
nodeID
!=
cj
->
nodeID
;
...
...
@@ -1341,6 +1424,14 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if
(
l
==
NULL
)
error
(
"Missing link to send_ti task."
);
scheduler_activate
(
s
,
l
->
t
);
}
}
else
if
(
t
->
type
==
task_type_pair
)
{
scheduler_activate
(
s
,
ci
->
drift
);
scheduler_activate
(
s
,
cj
->
drift
);
}
#else
if
(
t
->
type
==
task_type_pair
)
{
scheduler_activate
(
s
,
ci
->
drift
);
scheduler_activate
(
s
,
cj
->
drift
);
}
#endif
}
...
...
@@ -1355,7 +1446,6 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate
(
s
,
l
->
t
);
if
(
c
->
extra_ghost
!=
NULL
)
scheduler_activate
(
s
,
c
->
extra_ghost
);
if
(
c
->
ghost
!=
NULL
)
scheduler_activate
(
s
,
c
->
ghost
);
if
(
c
->
init
!=
NULL
)
scheduler_activate
(
s
,
c
->
init
);
if
(
c
->
init_grav
!=
NULL
)
scheduler_activate
(
s
,
c
->
init_grav
);
if
(
c
->
drift
!=
NULL
)
scheduler_activate
(
s
,
c
->
drift
);
if
(
c
->
kick1
!=
NULL
)
scheduler_activate
(
s
,
c
->
kick1
);
...
...
@@ -1409,7 +1499,9 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
/* Drift from the last time the cell was drifted to the current time */
const
double
dt
=
(
ti_current
-
ti_old
)
*
timeBase
;
float
dx_max
=
0
.
f
,
dx2_max
=
0
.
f
,
cell_h_max
=
0
.
f
;
float
dx_max
=
0
.
f
,
dx2_max
=
0
.
f
;
float
dx_max_sort
=
0
.
0
f
,
dx2_max_sort
=
0
.
f
;
float
cell_h_max
=
0
.
f
;
/* Check that we are actually going to move forward. */
if
(
ti_current
<
ti_old
)
error
(
"Attempt to drift to the past"
);
...
...
@@ -1421,8 +1513,13 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
for
(
int
k
=
0
;
k
<
8
;
k
++
)
if
(
c
->
progeny
[
k
]
!=
NULL
)
{
struct
cell
*
cp
=
c
->
progeny
[
k
];
/* Collect */
cell_drift_particles
(
cp
,
e
);
/* Update */
dx_max
=
max
(
dx_max
,
cp
->
dx_max
);
dx_max_sort
=
max
(
dx_max_sort
,
cp
->
dx_max_sort
);
cell_h_max
=
max
(
cell_h_max
,
cp
->
h_max
);
}
...
...
@@ -1443,6 +1540,11 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
gp
->
x_diff
[
1
]
*
gp
->
x_diff
[
1
]
+
gp
->
x_diff
[
2
]
*
gp
->
x_diff
[
2
];
dx2_max
=
max
(
dx2_max
,
dx2
);
/* Init gravity force fields. */
if
(
gpart_is_active
(
gp
,
e
))
{
gravity_init_gpart
(
gp
);
}
}
/* Loop over all the gas particles in the cell */
...
...
@@ -1464,9 +1566,18 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
xp
->
x_diff
[
1
]
*
xp
->
x_diff
[
1
]
+
xp
->
x_diff
[
2
]
*
xp
->
x_diff
[
2
];
dx2_max
=
max
(
dx2_max
,
dx2
);
const
float
dx2_sort
=
xp
->
x_diff_sort
[
0
]
*
xp
->
x_diff_sort
[
0
]
+
xp
->
x_diff_sort
[
1
]
*
xp
->
x_diff_sort
[
1
]
+
xp
->
x_diff_sort
[
2
]
*
xp
->
x_diff_sort
[
2
];
dx2_max_sort
=
max
(
dx2_max_sort
,
dx2_sort
);
/* Maximal smoothing length */
cell_h_max
=
max
(
cell_h_max
,
p
->
h
);
/* Get ready for a density calculation */
if
(
part_is_active
(
p
,
e
))
{
hydro_init_part
(
p
,
&
e
->
s
->
hs
);
}
}
/* Loop over all the star particles in the cell */
...
...
@@ -1484,16 +1595,19 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
/* Now, get the maximal particle motion from its square */
dx_max
=
sqrtf
(
dx2_max
);
dx_max_sort
=
sqrtf
(
dx2_max_sort
);
}
else
{
cell_h_max
=
c
->
h_max
;
dx_max
=
c
->
dx_max
;
dx_max_sort
=
c
->
dx_max_sort
;
}
/* Store the values */
c
->
h_max
=
cell_h_max
;
c
->
dx_max
=
dx_max
;
c
->
dx_max_sort
=
dx_max_sort
;
/* Update the time of the last drift */
c
->
ti_old
=
ti_current
;
...
...
src/cell.h
View file @
0167dd62
...
...
@@ -148,9 +148,6 @@ struct cell {
/*! Linked list of the tasks computing this cell's gravity forces. */