From e001fa4d5c2588e87bca908b63e57bf732361d0e Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Tue, 1 Aug 2017 09:17:07 +0100
Subject: [PATCH] Added a parameter to the gravity properties to contain the
distance below which all computed forces are Newtonian
---
examples/EAGLE_100/eagle_100.yml | 2 --
examples/EAGLE_12/eagle_12.yml | 2 --
examples/EAGLE_25/eagle_25.yml | 2 --
examples/EAGLE_50/eagle_50.yml | 6 ++++++
examples/UniformDMBox/uniformBox.yml | 3 +--
examples/parameter_example.yml | 11 ++++++-----
src/gravity_properties.c | 19 ++++++++++++-------
src/gravity_properties.h | 11 +++++++++--
src/runner_doiact_grav.h | 4 +++-
src/tools.c | 2 +-
10 files changed, 38 insertions(+), 24 deletions(-)
diff --git a/examples/EAGLE_100/eagle_100.yml b/examples/EAGLE_100/eagle_100.yml
index 54e56794c8..1ea1518825 100644
--- a/examples/EAGLE_100/eagle_100.yml
+++ b/examples/EAGLE_100/eagle_100.yml
@@ -28,8 +28,6 @@ Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
epsilon: 0.0001 # Softening length (in internal units).
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- a_smooth: 1.25 # PM smoothing scale in top-level cell-size units
- r_cut: 4.5 # Distance beyong wwhich FMM forces are zero in top-level cell-size units
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml
index 4bb04e757f..fd9569c50f 100644
--- a/examples/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_12/eagle_12.yml
@@ -28,8 +28,6 @@ Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
epsilon: 0.0001 # Softening length (in internal units).
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- a_smooth: 1.25 # PM smoothing scale in top-level cell-size units
- r_cut: 4.5 # Distance beyong wwhich FMM forces are zero in top-level cell-size un
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml
index 8920fb89a9..5dee9dad0b 100644
--- a/examples/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_25/eagle_25.yml
@@ -28,8 +28,6 @@ Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
epsilon: 0.0001 # Softening length (in internal units).
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- a_smooth: 1.25 # PM smoothing scale in top-level cell-size units
- r_cut: 4.5 # Distance beyong wwhich FMM forces are zero in top-level cell-size units
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_50/eagle_50.yml b/examples/EAGLE_50/eagle_50.yml
index b84b1eb7c3..898c28935a 100644
--- a/examples/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_50/eagle_50.yml
@@ -23,6 +23,12 @@ Snapshots:
Statistics:
delta_time: 1e-2 # Time between statistics output
+# Parameters for the self-gravity scheme
+Gravity:
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ epsilon: 0.0001 # Softening length (in internal units).
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
+
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/UniformDMBox/uniformBox.yml b/examples/UniformDMBox/uniformBox.yml
index dd8ea53a20..e59d677b30 100644
--- a/examples/UniformDMBox/uniformBox.yml
+++ b/examples/UniformDMBox/uniformBox.yml
@@ -28,8 +28,7 @@ Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
epsilon: 0.00001 # Softening length (in internal units).
- a_smooth: 1.25 # PM smoothing scale in top-level cell-size units
- r_cut: 4.5 # Distance beyong wwhich FMM forces are zero in top-level cell-size un
+
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 7ca15576bc..9c3cee7630 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -51,11 +51,12 @@ SPH:
# Parameters for the self-gravity scheme
Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- theta: 0.7 # Opening angle (Multipole acceptance criterion)
- epsilon: 0.1 # Softening length (in internal units).
- a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
- r_cut: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
+ epsilon: 0.1 # Softening length (in internal units).
+ a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
+ r_cut_max: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
+ r_cut_min: 0.1 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
# Parameters related to the initial conditions
InitialConditions:
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index b1098888b9..18cf044434 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -33,7 +33,8 @@
#include "kernel_gravity.h"
#define gravity_props_default_a_smooth 1.25f
-#define gravity_props_default_r_cut 4.5f
+#define gravity_props_default_r_cut_max 4.5f
+#define gravity_props_default_r_cut_min 0.1f
void gravity_props_init(struct gravity_props *p,
const struct swift_params *params) {
@@ -41,8 +42,10 @@ void gravity_props_init(struct gravity_props *p,
/* Tree-PM parameters */
p->a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth",
gravity_props_default_a_smooth);
- p->r_cut = parser_get_opt_param_float(params, "Gravity:r_cut",
- gravity_props_default_r_cut);
+ p->r_cut_max = parser_get_opt_param_float(params, "Gravity:r_cut_max",
+ gravity_props_default_r_cut_max);
+ p->r_cut_min = parser_get_opt_param_float(params, "Gravity:r_cut_min",
+ gravity_props_default_r_cut_min);
/* Time integration */
p->eta = parser_get_param_float(params, "Gravity:eta");
@@ -69,9 +72,10 @@ void gravity_props_print(const struct gravity_props *p) {
message("Self-gravity softening: epsilon=%.4f (Plummer equivalent: %.4f)",
p->epsilon, p->epsilon / 3.);
- message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth);
+ message("Self-gravity mesh smoothing-scale: a_smooth=%f", p->a_smooth);
- message("Self-gravity MM cut-off: r_cut=%f", p->r_cut);
+ message("Self-gravity tree cut-off: r_cut_max=%f", p->r_cut_max);
+ message("Self-gravity truncation cut-off: r_cut_min=%f", p->r_cut_min);
}
#if defined(HAVE_HDF5)
@@ -84,7 +88,8 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
p->epsilon / 3.);
io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
- io_write_attribute_f(h_grpgrav, "MM a_smooth", p->a_smooth);
- io_write_attribute_f(h_grpgrav, "MM r_cut", p->r_cut);
+ io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
+ io_write_attribute_f(h_grpgrav, "Mesh r_cut_max", p->r_cut_max);
+ io_write_attribute_f(h_grpgrav, "Mesh r_cut_min", p->r_cut_min);
}
#endif
diff --git a/src/gravity_properties.h b/src/gravity_properties.h
index be26f0d1d2..2a5e4cb1e0 100644
--- a/src/gravity_properties.h
+++ b/src/gravity_properties.h
@@ -34,9 +34,16 @@
*/
struct gravity_props {
- /* Tree-PM parameters */
+ /*! Mesh smoothing scale in units of top-level cell size */
float a_smooth;
- float r_cut;
+
+ /*! Distance below which the truncated mesh force is Newtonian in units of
+ * a_smooth */
+ float r_cut_min;
+
+ /*! Distance above which the truncated mesh force is negligible in units of
+ * a_smooth */
+ float r_cut_max;
/*! Time integration dimensionless multiplier */
float eta;
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 1de491fa01..acf6478a09 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -664,7 +664,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
const double cell_width = s->width[0];
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const double theta_crit_inv = props->theta_crit_inv;
- const double max_distance = props->a_smooth * props->r_cut * cell_width;
+ const double max_distance = props->a_smooth * props->r_cut_max * cell_width;
const double max_distance2 = max_distance * max_distance;
struct gravity_tensors *const mi = ci->multipole;
const double CoM[3] = {mi->CoM[0], mi->CoM[1], mi->CoM[2]};
@@ -701,9 +701,11 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
const double dz = nearest(CoM[2] - mj->CoM[2], dim[2]);
const double r2 = dx * dx + dy * dy + dz * dz;
+ /* Are we beyond the distance where the truncated forces are 0 ?*/
if (r2 > max_distance2) {
#ifdef SWIFT_DEBUG_CHECKS
+ /* Need to account for the interactions we missed */
mi->pot.num_interacted += mj->m_pole.num_gpart;
#endif
continue;
diff --git a/src/tools.c b/src/tools.c
index d701afc5fd..0f423d5aad 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -755,7 +755,7 @@ void gravity_n2(struct gpart *gparts, const int gcount,
const struct gravity_props *gravity_properties, float rlr) {
const float rlr_inv = 1. / rlr;
- const float r_cut = gravity_properties->r_cut;
+ const float r_cut = gravity_properties->r_cut_max;
const float max_d = r_cut * rlr;
const float max_d2 = max_d * max_d;
--
GitLab