From 47bfee49c78b899545d37316d7f7ed3b62aa7e1c Mon Sep 17 00:00:00 2001
From: rttw52 <stuart.mcalpine@helsink.fi>
Date: Fri, 13 Mar 2020 23:34:46 +0000
Subject: [PATCH] Sorry I messed up the git so it's one big commit.
Updated RTD as per your recomendations.
Added ForceChecks to paramater_example.yml
Changed brute_force_gravity_flag to force_checks_snapshot_flag
Removed epsilon from gravity_exact_force_file_exits
Made functions in src/kernel_long_gravity.h double
---
.../source/GettingStarted/special_modes.rst | 11 +++---
.../ParameterFiles/parameter_description.rst | 19 +++++-----
examples/parameter_example.yml | 6 +++-
src/engine.c | 36 ++++++++++---------
src/engine.h | 4 +--
src/gravity.c | 7 +---
src/kernel_long_gravity.h | 26 +++++++-------
7 files changed, 57 insertions(+), 52 deletions(-)
diff --git a/doc/RTD/source/GettingStarted/special_modes.rst b/doc/RTD/source/GettingStarted/special_modes.rst
index 47d6cbd5c5..ce89d857ab 100644
--- a/doc/RTD/source/GettingStarted/special_modes.rst
+++ b/doc/RTD/source/GettingStarted/special_modes.rst
@@ -53,8 +53,8 @@ simply the Newtonian sum, i.e.,
To run SWIFT in this mode, configure the code with
``--enable-gravity-force-checks=N``, which means that the exact forces will be
-computed for every :math:`N^{th}` particle (i.e., to compute the exact forces
-for all particles set ``N=1``).
+computed for every :math:`N^{th}` particle based on their ID (i.e., to compute
+the exact forces for all particles set ``N=1``).
Two `.dat` files will be output during each timestep, one containing the forces
(really it is the accelerations that are stored) as computed by ``_swift_``, and
@@ -66,5 +66,8 @@ the number of particles contributing to each force component is also stored
(these numbers will add up to :math:`n-1`).
This mode will slow down the code *considerably*, and it is not recommended to
-run on problems with more than :math:`10^{5}` particles. This mode must be run
-on a single node/rank and is only designed for pure gravity tests (i.e., DMO).
+be run on problems with more than :math:`10^{5}` particles when
+``--enable-gravity-force-checks=1``. For larger runs, sampling a sub-set of
+particles via the argument ``N`` of the configuration option is recommended.
+This mode must be run on a single node/rank, and is primarily designed for pure
+gravity tests (i.e., DMO).
diff --git a/doc/RTD/source/ParameterFiles/parameter_description.rst b/doc/RTD/source/ParameterFiles/parameter_description.rst
index ab0c1cc349..1fade92ba6 100644
--- a/doc/RTD/source/ParameterFiles/parameter_description.rst
+++ b/doc/RTD/source/ParameterFiles/parameter_description.rst
@@ -1258,19 +1258,20 @@ Showing all the parameters for a basic cosmologica test-case, one would have:
Gravity Force Checks
--------------------
-By default, when configured with ``--enable-gravity-force-checks`` the "exact"
-forces will be computed for all active particles during each timestep.
+By default, when the code is configured with ``--enable-gravity-force-checks``,
+the "exact" forces of all active gparts are computed during each timestep.
-To give a bit more control over this, you can select to only compute the exact
-forces during the timesteps that all the particles are active, and/or only at
-the timesteps when a snapshot is being dumped, i.e.,
+To give a bit more control over this, you can select to only perform the exact
+force computation during the timesteps that all gparts are active, and/or only
+at the timesteps when a snapshot is being dumped, i.e.,
.. code:: YAML
ForceChecks:
- only_when_all_active: 1 # Only compute exact forces during timesteps when all particles are active.
+ only_when_all_active: 1 # Only compute exact forces during timesteps when all gparts are active.
only_at_snapshots: 1 # Only compute exact forces during timesteps when a snapshot is being dumped.
-Note that if ``only_when_all_active:1`` and ``only_when_all_active:1``, the
-exact forces will be computed during the first timestep after the snapshot
-output time when all particles are active.
+If ``only_when_all_active:1`` and ``only_when_all_active:1`` are enabled
+together, and all the gparts are not active during the timestep of the snapshot
+dump, the exact forces computation is performed on the first timestep at which
+all the gparts are active after that snapshot output timestep.
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index cabc94d2bb..6a3e1fa9ee 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -77,7 +77,11 @@ Gravity:
dithering: 1 # (Optional) Activate the dithering of the gravity mesh at every rebuild (this is the default value).
dithering_ratio: 1.0 # (Optional) Magnitude of each component of the dithering vector in units of the top-level cell sizes (this is the default value).
-# Parameters for the Friends-Of-Friends algorithm
+# Parameters when running with SWIFT_GRAVITY_FORCE_CHECKS
+ForceChecks:
+ only_when_all_active: 1 # (Optional) Only compute exact forces during timesteps when all gparts are active (default: 0).
+ only_at_snapshots: 1 # (Optional) Only compute exact forces during timesteps when a snapshot is being dumped (default: 0).
+
FOF:
basename: fof_output # Filename for the FOF outputs (Unused when FoF is only run to seed BHs).
scale_factor_first: 0.91 # Scale-factor of first FoF black hole seeding calls (needed for cosmological runs).
diff --git a/src/engine.c b/src/engine.c
index d0f367cd73..ae826518b2 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2165,7 +2165,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
e->wallclock_time = (float)clocks_diff(&time1, &time2);
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
- e->brute_force_gravity_flag = 0;
+ e->force_checks_snapshot_flag = 0;
#endif
if (e->verbose) message("took %.3f %s.", e->wallclock_time, clocks_getunit());
@@ -2362,12 +2362,13 @@ void engine_step(struct engine *e) {
if ((e->all_gparts_active && e->force_checks_only_all_active) ||
!e->force_checks_only_all_active) {
/* Is this a snapshot timestep (and the option is selected)? */
- if ((e->brute_force_gravity_flag == 1 &&
- e->force_checks_only_at_snapshots) ||
- !e->force_checks_only_at_snapshots) {
-
- /* Do checks */
- gravity_exact_force_compute(e->s, e);
+ if ((e->force_checks_snapshot_flag &&
+ e->force_checks_only_at_snapshots) ||
+ !e->force_checks_only_at_snapshots) {
+
+ /* Do checks */
+ gravity_exact_force_compute(e->s, e);
+ }
}
}
}
@@ -2395,15 +2396,16 @@ void engine_step(struct engine *e) {
if ((e->all_gparts_active && e->force_checks_only_all_active) ||
!e->force_checks_only_all_active) {
/* Is this a snapshot timestep (and the option is selected)? */
- if ((e->brute_force_gravity_flag == 1 &&
- e->force_checks_only_at_snapshots) ||
- !e->force_checks_only_at_snapshots) {
-
- /* Do checks */
- gravity_exact_force_check(e->s, e, 1e-1);
-
- /* Reset flag waiting for next output time */
- e->brute_force_gravity_flag = 0;
+ if ((e->force_checks_snapshot_flag &&
+ e->force_checks_only_at_snapshots) ||
+ !e->force_checks_only_at_snapshots) {
+
+ /* Do checks */
+ gravity_exact_force_check(e->s, e, 1e-1);
+
+ /* Reset flag waiting for next output time */
+ e->force_checks_snapshot_flag = 0;
+ }
}
}
}
@@ -2535,7 +2537,7 @@ void engine_check_for_dumps(struct engine *e) {
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
/* Indicate we are allowed to do a brute force calculation now */
- e->brute_force_gravity_flag = 1;
+ e->force_checks_snapshot_flag = 1;
#endif
/* Do we want a corresponding VELOCIraptor output? */
diff --git a/src/engine.h b/src/engine.h
index d34b9ea1b5..3e35dbb0dd 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -495,8 +495,8 @@ struct engine {
/* Are all gparts active this timestep? */
int all_gparts_active;
- /* Are we doing brute force calculations this timestep? */
- int brute_force_gravity_flag;
+ /* Flag to tell brute force checks a snapshot was recently written. */
+ int force_checks_snapshot_flag;
#endif
};
diff --git a/src/gravity.c b/src/gravity.c
index 5e0a92503a..7e9c4d02f0 100644
--- a/src/gravity.c
+++ b/src/gravity.c
@@ -430,7 +430,7 @@ int gravity_exact_force_file_exits(const struct engine *e) {
char line[100];
char dummy1[10], dummy2[10];
- double epsilon, newton_G;
+ double newton_G;
int N, periodic;
/* Reads file header */
if (fgets(line, 100, file) != line) error("Problem reading title");
@@ -438,17 +438,12 @@ int gravity_exact_force_file_exits(const struct engine *e) {
sscanf(line, "%s %s %le", dummy1, dummy2, &newton_G);
if (fgets(line, 100, file) != line) error("Problem reading N");
sscanf(line, "%s %s %d", dummy1, dummy2, &N);
- if (fgets(line, 100, file) != line) error("Problem reading epsilon");
- sscanf(line, "%s %s %le", dummy1, dummy2, &epsilon);
if (fgets(line, 100, file) != line) error("Problem reading BC");
sscanf(line, "%s %s %d", dummy1, dummy2, &periodic);
fclose(file);
/* Check whether it matches the current parameters */
if (N == SWIFT_GRAVITY_FORCE_CHECKS && periodic == e->s->periodic &&
- (fabs(epsilon - gravity_get_softening(0, e->gravity_properties)) /
- epsilon <
- 1e-5) &&
(fabs(newton_G - e->physical_constants->const_newton_G) / newton_G <
1e-5)) {
return 1;
diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h
index 25d1d64338..ae8a828cae 100644
--- a/src/kernel_long_gravity.h
+++ b/src/kernel_long_gravity.h
@@ -222,25 +222,25 @@ kernel_long_grav_force_eval_double(const double u, double *const W) {
const double one_over_sqrt_pi = ((double)(M_2_SQRTPI * 0.5));
- const double arg1 = u * 0.5f;
- const double arg2 = -arg1 * arg1;
+ const double arg1 = u * 0.5;
+ const double arg2 = -arg1 * arg1;
- const double term1 = approx_erfcf(arg1);
- const double term2 = u * one_over_sqrt_pi * expf(arg2);
+ const double term1 = erfc(arg1);
+ const double term2 = u * one_over_sqrt_pi * exp(arg2);
*W = term1 + term2;
#else
- const double x = 2.f * u;
- const double exp_x = expf(x); // good_approx_expf(x);
- const double alpha = 1.f / (1.f + exp_x);
+ const double x = 2. * u;
+ const double exp_x = exp(x);
+ const double alpha = 1. / (1. + exp_x);
- /* We want 2*(x*alpha - x*alpha^2 - exp(x)*alpha + 1) */
- *W = 1.f - alpha;
- *W = *W * x - exp_x;
- *W = *W * alpha + 1.f;
- *W *= 2.f;
-#endif
+ /* We want 2*(x*alpha - x*alpha^2 - exp(x)*alpha + 1) */
+ *W = 1. - alpha;
+ *W = *W * x - exp_x;
+ *W = *W * alpha + 1.;
+ *W *= 2.;
+ #endif
#endif
}
--
GitLab