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
090a5955
Commit
090a5955
authored
Aug 15, 2017
by
Matthieu Schaller
Browse files
Initial implementation of mass-dependent multipole acceptance criterion
parent
e0a98486
Changes
7
Hide whitespace changes
Inline
Side-by-side
examples/EAGLE_6/eagle_6.yml
View file @
090a5955
...
...
@@ -30,7 +30,7 @@ Statistics:
Gravity
:
eta
:
0.025
# Constant dimensionless multiplier for time integration.
theta
:
0.7
# Opening angle (Multipole acceptance criterion)
epsilon
:
0.00
0
1
# Softening length (in internal units).
epsilon
:
0.001
# Softening length (in internal units).
# Parameters for the hydrodynamics scheme
SPH
:
...
...
src/cell.c
View file @
090a5955
...
...
@@ -1679,7 +1679,7 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
const
struct
engine
*
e
=
sp
->
e
;
const
int
periodic
=
sp
->
periodic
;
const
double
dim
[
3
]
=
{
sp
->
dim
[
0
],
sp
->
dim
[
1
],
sp
->
dim
[
2
]};
const
double
theta_crit
=
e
->
gravity_properties
->
theta_crit
;
const
double
theta_crit
2
=
e
->
gravity_properties
->
theta_crit
2
;
/* Self interaction? */
if
(
cj
==
NULL
)
{
...
...
@@ -1733,7 +1733,7 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
const
double
r2
=
dx
*
dx
+
dy
*
dy
+
dz
*
dz
;
/* Can we use multipoles ? */
if
(
gravity_multipole_accept
(
multi_i
,
multi_j
,
theta_crit
,
r2
))
{
if
(
gravity_multipole_accept
(
multi_i
,
multi_j
,
theta_crit
2
,
r2
))
{
/* Ok, no need to drift anything */
return
;
...
...
src/engine.c
View file @
090a5955
...
...
@@ -1726,7 +1726,7 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
const
int
cdim
[
3
]
=
{
s
->
cdim
[
0
],
s
->
cdim
[
1
],
s
->
cdim
[
2
]};
const
int
cdim_ghost
[
3
]
=
{
s
->
cdim
[
0
]
/
4
+
1
,
s
->
cdim
[
1
]
/
4
+
1
,
s
->
cdim
[
2
]
/
4
+
1
};
const
double
theta_crit
=
e
->
gravity_properties
->
theta_crit
;
const
double
theta_crit
2
=
e
->
gravity_properties
->
theta_crit
2
;
struct
cell
*
cells
=
s
->
cells_top
;
const
int
n_ghosts
=
cdim_ghost
[
0
]
*
cdim_ghost
[
1
]
*
cdim_ghost
[
2
]
*
2
;
...
...
@@ -1798,7 +1798,7 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
const
double
r2
=
dx
*
dx
+
dy
*
dy
+
dz
*
dz
;
/* Are the cells too close for a MM interaction ? */
if
(
!
gravity_multipole_accept_rebuild
(
multi_i
,
multi_j
,
theta_crit
,
if
(
!
gravity_multipole_accept_rebuild
(
multi_i
,
multi_j
,
theta_crit
2
,
r2
))
{
/* Ok, we need to add a direct pair calculation */
...
...
src/gravity_properties.c
View file @
090a5955
...
...
@@ -52,6 +52,8 @@ void gravity_props_init(struct gravity_props *p,
/* Opening angle */
p
->
theta_crit
=
parser_get_param_double
(
params
,
"Gravity:theta"
);
if
(
p
->
theta_crit
>=
1
.)
error
(
"Theta too large. FMM won't converge."
);
p
->
theta_crit2
=
p
->
theta_crit
*
p
->
theta_crit
;
p
->
theta_crit_inv
=
1
.
/
p
->
theta_crit
;
/* Softening lengths */
...
...
src/gravity_properties.h
View file @
090a5955
...
...
@@ -51,6 +51,9 @@ struct gravity_props {
/*! Tree opening angle (Multipole acceptance criterion) */
double
theta_crit
;
/*! Square of opening angle */
double
theta_crit2
;
/*! Inverse of opening angle */
double
theta_crit_inv
;
...
...
src/multipole.h
View file @
090a5955
...
...
@@ -234,6 +234,7 @@ INLINE static void gravity_drift(struct gravity_tensors *m, double dt) {
* @brief Zeroes all the fields of a field tensor
*
* @param l The field tensor.
* @param ti_current The current (integer) time (for debugging only).
*/
INLINE
static
void
gravity_field_tensors_init
(
struct
grav_tensor
*
l
,
integertime_t
ti_current
)
{
...
...
@@ -246,7 +247,7 @@ INLINE static void gravity_field_tensors_init(struct grav_tensor *l,
}
/**
* @brief Adds field tensr
s
to other one
s
(i.e. does la += lb).
* @brief Adds
a
field tens
o
r to
an
other one (i.e. does la += lb).
*
* @param la The gravity tensors to add to.
* @param lb The gravity tensors to add.
...
...
@@ -2174,6 +2175,9 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
l_b
->
F_102
+=
m_a
->
M_000
*
D_soft_102
(
dx
,
dy
,
dz
,
r
,
eps_inv
);
l_b
->
F_012
+=
m_a
->
M_000
*
D_soft_012
(
dx
,
dy
,
dz
,
r
,
eps_inv
);
l_b
->
F_111
+=
m_a
->
M_000
*
D_soft_111
(
dx
,
dy
,
dz
,
r
,
eps_inv
);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif
}
}
...
...
@@ -2654,22 +2658,24 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
*
* @param ma The #multipole of the first #cell.
* @param mb The #multipole of the second #cell.
* @param theta_crit The critical opening angle.
* @param theta_crit
2
The
square of the
critical opening angle.
* @param r2 Square of the distance (periodically wrapped) between the
* multipoles.
*/
__attribute__
((
always_inline
))
INLINE
static
int
gravity_multipole_accept_rebuild
(
const
struct
gravity_tensors
*
const
ma
,
const
struct
gravity_tensors
*
const
mb
,
double
theta_crit
,
double
r2
)
{
double
theta_crit
2
,
double
r2
)
{
const
double
r_crit_a
=
ma
->
r_max_rebuild
*
theta_crit
;
const
double
r_crit_b
=
mb
->
r_max_rebuild
*
theta_crit
;
const
double
r_crit_a
=
ma
->
r_max_rebuild
;
const
double
r_crit_b
=
mb
->
r_max_rebuild
;
const
double
size
=
r_crit_a
+
r_crit_b
;
const
double
size2
=
size
*
size
;
// MATTHIEU: Make this mass-dependent ?
/* Multipole acceptance criterion (Dehnen 2002, eq.10) */
return
(
r2
>
(
r_crit_a
+
r_crit_b
)
*
(
r_crit_a
+
r_crit_b
)
);
return
(
r2
*
theta_crit2
>
size2
);
}
/**
...
...
@@ -2681,21 +2687,23 @@ gravity_multipole_accept_rebuild(const struct gravity_tensors *const ma,
*
* @param ma The #multipole of the first #cell.
* @param mb The #multipole of the second #cell.
* @param theta_crit The critical opening angle.
* @param theta_crit
2
The
square of the
critical opening angle.
* @param r2 Square of the distance (periodically wrapped) between the
* multipoles.
*/
__attribute__
((
always_inline
))
INLINE
static
int
gravity_multipole_accept
(
const
struct
gravity_tensors
*
const
ma
,
const
struct
gravity_tensors
*
const
mb
,
double
theta_crit
,
double
r2
)
{
const
struct
gravity_tensors
*
const
mb
,
double
theta_crit
2
,
double
r2
)
{
const
double
r_crit_a
=
ma
->
r_max
*
theta_crit
;
const
double
r_crit_b
=
mb
->
r_max
*
theta_crit
;
const
double
r_crit_a
=
ma
->
r_max
;
const
double
r_crit_b
=
mb
->
r_max
;
const
double
size
=
r_crit_a
+
r_crit_b
;
const
double
size2
=
size
*
size
;
// MATTHIEU: Make this mass-dependent ?
/* Multipole acceptance criterion (Dehnen 2002, eq.10) */
return
(
r2
>
(
r_crit_a
+
r_crit_b
)
*
(
r_crit_a
+
r_crit_b
)
);
return
(
r2
*
theta_crit2
>
size2
);
}
#endif
/* SWIFT_MULTIPOLE_H */
src/runner_doiact_grav.h
View file @
090a5955
...
...
@@ -1412,7 +1412,7 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
const
double
cell_width
=
s
->
width
[
0
];
const
double
dim
[
3
]
=
{
s
->
dim
[
0
],
s
->
dim
[
1
],
s
->
dim
[
2
]};
const
struct
gravity_props
*
props
=
e
->
gravity_properties
;
const
double
theta_crit
=
props
->
theta_crit
;
const
double
theta_crit
2
=
props
->
theta_crit
2
;
const
double
max_distance
=
props
->
a_smooth
*
props
->
r_cut_max
*
cell_width
;
const
double
max_distance2
=
max_distance
*
max_distance
;
...
...
@@ -1473,7 +1473,7 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
* option... */
/* Can we use M-M interactions ? */
if
(
gravity_multipole_accept
(
multi_i
,
multi_j
,
theta_crit
,
r2
))
{
if
(
gravity_multipole_accept
(
multi_i
,
multi_j
,
theta_crit
2
,
r2
))
{
/* MATTHIEU: make a symmetric M-M interaction function ! */
runner_dopair_grav_mm
(
r
,
ci
,
cj
);
...
...
@@ -1638,7 +1638,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
const
int
periodic
=
s
->
periodic
;
const
double
cell_width
=
s
->
width
[
0
];
const
double
dim
[
3
]
=
{
s
->
dim
[
0
],
s
->
dim
[
1
],
s
->
dim
[
2
]};
const
double
theta_crit
=
props
->
theta_crit
;
const
double
theta_crit
2
=
props
->
theta_crit
2
;
const
double
max_distance
=
props
->
a_smooth
*
props
->
r_cut_max
*
cell_width
;
const
double
max_distance2
=
max_distance
*
max_distance
;
...
...
@@ -1697,7 +1697,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
}
/* Check the multipole acceptance criterion */
if
(
gravity_multipole_accept
(
multi_i
,
multi_j
,
theta_crit
,
r2
))
{
if
(
gravity_multipole_accept
(
multi_i
,
multi_j
,
theta_crit
2
,
r2
))
{
/* Go for a (non-symmetric) M-M calculation */
runner_dopair_grav_mm
(
r
,
ci
,
cj
);
...
...
@@ -1720,7 +1720,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
const
double
r2_rebuild
=
dx
*
dx
+
dy
*
dy
+
dz
*
dz
;
/* Is the criterion violated now but was OK at the last rebuild ? */
if
(
gravity_multipole_accept_rebuild
(
multi_i
,
multi_j
,
theta_crit
,
if
(
gravity_multipole_accept_rebuild
(
multi_i
,
multi_j
,
theta_crit
2
,
r2_rebuild
))
{
/* Alright, we have to take charge of that pair in a different way. */
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment