From e8c8f69a16de9c55ff6885e12e463536a5ba808c Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sun, 1 May 2016 15:18:46 +0200
Subject: [PATCH] SPH properties are now read from the parameter file
---
examples/UniformBox/uniformBox.yml | 9 +++--
examples/main.c | 6 ++-
src/Makefile.am | 4 +-
src/cell.c | 1 +
src/const.h | 24 +++++-------
src/engine.c | 16 ++++----
src/engine.h | 5 +++
src/hydro/Gadget2/hydro.h | 16 +++-----
src/hydro/Gadget2/hydro_io.h | 18 ++++-----
src/hydro_properties.c | 61 ++++++++++++++++++++++++++++++
src/hydro_properties.h | 56 +++++++++++++++++++++++++++
src/kernel_hydro.h | 8 +---
src/potentials.h | 1 +
src/runner.c | 22 ++++++++---
src/swift.h | 1 +
15 files changed, 187 insertions(+), 61 deletions(-)
create mode 100644 src/hydro_properties.c
create mode 100644 src/hydro_properties.h
diff --git a/examples/UniformBox/uniformBox.yml b/examples/UniformBox/uniformBox.yml
index 5ff07586f3..9b7139b7c4 100644
--- a/examples/UniformBox/uniformBox.yml
+++ b/examples/UniformBox/uniformBox.yml
@@ -34,12 +34,13 @@ Snapshots:
# Parameters for the hydrodynamics scheme
SPH:
- resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
- delta_neighbours: 1. # The tolerance for the targetted number of neighbours.
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
max_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
max_smoothing_length: 0.1 # Maximal smoothing length allowed (in internal units).
-
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ max_volume_change: 2. # Maximal allowed change of volume over one time-step
+
# Parameters related to the initial conditions
InitialConditions:
file_name: ./uniformBox.hdf5 # The file to read
diff --git a/examples/main.c b/examples/main.c
index f4dd3d93e3..e5417c6286 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -302,6 +302,10 @@ int main(int argc, char *argv[]) {
phys_const_print(&prog_const);
}
+ /* Initialise the hydro properties */
+ struct hydro_props hydro_properties;
+ hydro_props_init(&hydro_properties, params);
+
/* Initialise the external potential properties */
struct external_potential potential;
if (with_external_gravity) potential_init(params, &us, &potential);
@@ -414,7 +418,7 @@ int main(int argc, char *argv[]) {
if (myrank == 0) clocks_gettime(&tic);
struct engine e;
engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, engine_policies,
- talking, &prog_const, &potential);
+ talking, &prog_const, &hydro_properties, &potential);
if (myrank == 0) {
clocks_gettime(&toc);
message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
diff --git a/src/Makefile.am b/src/Makefile.am
index 57ca1e5974..f50f5574d1 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -36,14 +36,14 @@ endif
include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h \
- physical_constants.h physical_constants_cgs.h potentials.h version.h
+ physical_constants.h physical_constants_cgs.h potentials.h version.h hydro_properties.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
units.c common_io.c single_io.c multipole.c version.c map.c \
kernel_hydro.c kernel_gravity.c tools.c part.c partition.c clocks.c parser.c \
- physical_constants.c potentials.c
+ physical_constants.c potentials.c hydro_properties.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
diff --git a/src/cell.c b/src/cell.c
index 793a2e10c2..29072173b0 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -50,6 +50,7 @@
#include "atomic.h"
#include "error.h"
#include "gravity.h"
+#include "hydro_properties.h"
#include "hydro.h"
#include "space.h"
#include "timers.h"
diff --git a/src/const.h b/src/const.h
index dd3e34016c..2d84915fed 100644
--- a/src/const.h
+++ b/src/const.h
@@ -24,8 +24,7 @@
#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 0.8f
#define const_viscosity_alpha_min \
0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
#define const_viscosity_alpha_max \
@@ -38,19 +37,8 @@
1.f /* Value taken from (Price,2008), not used in legacy gadget mode */
/* Time integration constants. */
-#define const_cfl 0.1f
-#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 0.1f
-#define const_smoothing_max_iter 30
-#define CUBIC_SPLINE_KERNEL
-
/* Gravity stuff. */
#define const_theta_max \
0.57735f /* Opening criteria, which is the ratio of the \
@@ -65,12 +53,20 @@
#define const_iepsilon5 (const_iepsilon3* const_iepsilon2)
#define const_iepsilon6 (const_iepsilon3* const_iepsilon3)
+/* Kernel function to use */
+#define CUBIC_SPLINE_KERNEL
+//#define QUARTIC_SPLINE_KERNEL
+//#define QUINTIC_SPLINE_KERNEL
+//#define WENDLAND_C2_KERNEL
+//#define WENDLAND_C4_KERNEL
+//#define WENDLAND_C6_KERNEL
+
/* SPH variant to use */
//#define MINIMAL_SPH
#define GADGET2_SPH
//#define DEFAULT_SPH
/* Gravity properties */
-#define EXTERNAL_POTENTIAL_POINTMASS
+//#define EXTERNAL_POTENTIAL_POINTMASS
#endif /* SWIFT_CONST_H */
diff --git a/src/engine.c b/src/engine.c
index eabe688f29..af0edb9ff9 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2356,9 +2356,8 @@ void engine_dump_snapshot(struct engine *e) {
struct clocks_time time1, time2;
clocks_gettime(&time1);
- if (e->verbose)
- message("writing snapshot at t=%f", e->time);
-
+ if (e->verbose) message("writing snapshot at t=%f", e->time);
+
/* Dump... */
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
@@ -2409,6 +2408,7 @@ static bool hyperthreads_present(void) {
* @param policy The queuing policy to use.
* @param verbose Is this #engine talkative ?
* @param physical_constants The #phys_const used for this run.
+ * @param hydro The #hydro_props used for this run.
* @param potential The properties of the external potential.
*/
@@ -2416,6 +2416,7 @@ void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
int nr_threads, int policy, int verbose,
const struct phys_const *physical_constants,
+ const struct hydro_props *hydro,
const struct external_potential *potential) {
/* Clean-up everything */
@@ -2457,6 +2458,7 @@ void engine_init(struct engine *e, struct space *s,
e->count_step = 0;
e->wallclock_time = 0.f;
e->physical_constants = physical_constants;
+ e->hydro_properties = hydro;
e->external_potential = potential;
engine_rank = nodeID;
@@ -2563,12 +2565,8 @@ void engine_init(struct engine *e, struct space *s,
engine_print_policy(e);
/* Print information about the hydro scheme */
- if (e->policy & engine_policy_hydro) {
- if (e->nodeID == 0) message("Hydrodynamic scheme: %s.", SPH_IMPLEMENTATION);
- if (e->nodeID == 0)
- message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours.",
- kernel_name, kernel_nwneigh, const_delta_nwneigh);
- }
+ if (e->policy & engine_policy_hydro)
+ if (e->nodeID == 0) hydro_props_print(e->hydro_properties);
/* Check we have sensible time bounds */
if (e->timeBegin >= e->timeEnd)
diff --git a/src/engine.h b/src/engine.h
index f4bc3be00c..9ef7d57599 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -37,6 +37,7 @@
#include <stdio.h>
/* Includes. */
+#include "hydro_properties.h"
#include "lock.h"
#include "proxy.h"
#include "runner.h"
@@ -183,6 +184,9 @@ struct engine {
/* Physical constants definition */
const struct phys_const *physical_constants;
+ /* Properties of the hydro scheme */
+ const struct hydro_props *hydro_properties;
+
/* Properties of external gravitational potential */
const struct external_potential *external_potential;
};
@@ -195,6 +199,7 @@ void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
int nr_threads, int policy, int verbose,
const struct phys_const *physical_constants,
+ const struct hydro_props *hydro,
const struct external_potential *potential);
void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
unsigned int submask);
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 22c5734ed5..e84db11ad1 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -25,20 +25,16 @@
*
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
- struct part* p, struct xpart* xp) {
+ struct part* p, struct xpart* xp,
+ const struct hydro_props* hydro_properties) {
- /* Acceleration */
- float ac =
- sqrtf(p->a_hydro[0] * p->a_hydro[0] + p->a_hydro[1] * p->a_hydro[1] +
- p->a_hydro[2] * p->a_hydro[2]);
- ac = fmaxf(ac, 1e-30);
-
- const float dt_accel = sqrtf(2.f); // MATTHIEU
+ const float CFL_condition = hydro_properties->CFL_condition;
/* CFL condition */
- const float dt_cfl = 2.f * const_cfl * kernel_gamma * p->h / p->force.v_sig;
+ const float dt_cfl =
+ 2.f * kernel_gamma * CFL_condition * p->h / p->force.v_sig;
- return fminf(dt_cfl, dt_accel);
+ return dt_cfl;
}
/**
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index c1c59dfa49..b1db6a8938 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -104,10 +104,10 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Kernel function description */
writeAttribute_s(h_grpsph, "Kernel", kernel_name);
- writeAttribute_f(h_grpsph, "Kernel eta", const_eta_kernel);
- writeAttribute_f(h_grpsph, "Weighted N_ngb", kernel_nwneigh);
- writeAttribute_f(h_grpsph, "Delta N_ngb", const_delta_nwneigh);
- writeAttribute_f(h_grpsph, "Hydro gamma", const_hydro_gamma);
+ //writeAttribute_f(h_grpsph, "Kernel eta", const_eta_kernel);
+ //writeAttribute_f(h_grpsph, "Weighted N_ngb", kernel_nwneigh);
+ //writeAttribute_f(h_grpsph, "Delta N_ngb", const_delta_nwneigh);
+ //writeAttribute_f(h_grpsph, "Hydro gamma", const_hydro_gamma);
/* Viscosity and thermal conduction */
writeAttribute_s(h_grpsph, "Thermal Conductivity Model",
@@ -118,9 +118,9 @@ void writeSPHflavour(hid_t h_grpsph) {
writeAttribute_f(h_grpsph, "Viscosity beta", 3.f);
/* Time integration properties */
- writeAttribute_f(h_grpsph, "CFL parameter", const_cfl);
- writeAttribute_f(h_grpsph, "Maximal ln(Delta h) change over dt",
- const_ln_max_h_change);
- writeAttribute_f(h_grpsph, "Maximal Delta h change over dt",
- exp(const_ln_max_h_change));
+ //writeAttribute_f(h_grpsph, "CFL parameter", const_cfl);
+ //writeAttribute_f(h_grpsph, "Maximal ln(Delta h) change over dt",
+ // const_ln_max_h_change);
+ //writeAttribute_f(h_grpsph, "Maximal Delta h change over dt",
+ // exp(const_ln_max_h_change));
}
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
new file mode 100644
index 0000000000..03cdc41beb
--- /dev/null
+++ b/src/hydro_properties.c
@@ -0,0 +1,61 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 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/>.
+ *
+ ******************************************************************************/
+
+/* This object's header. */
+#include "hydro_properties.h"
+
+/* Standard headers */
+#include <math.h>
+
+/* Local headers. */
+#include "error.h"
+#include "kernel_hydro.h"
+#include "hydro.h"
+
+void hydro_props_init(struct hydro_props *p,
+ const struct swift_params *params) {
+
+ /* Kernel properties */
+ p->eta_neighbours = parser_get_param_float(params, "SPH:resolution_eta");
+ const float eta3 = p->eta_neighbours * p->eta_neighbours * p->eta_neighbours;
+ p->target_neighbours = 4.0 * M_PI * kernel_gamma3 * eta3 / 3.0;
+ p->delta_neighbours = parser_get_param_float(params, "SPH:delta_neighbours");
+
+ /* Ghost stuff */
+ p->max_smoothing_iterations =
+ parser_get_param_int(params, "SPH:max_ghost_iterations");
+
+ /* Time integration properties */
+ p->CFL_condition = parser_get_param_float(params, "SPH:CFL_condition");
+ const float max_volume_change =
+ parser_get_param_float(params, "SPH:max_volume_change");
+ p->log_max_h_change = log(powf(max_volume_change, 0.33333333333f));
+}
+
+void hydro_props_print(const struct hydro_props *p) {
+
+ message("Hydrodynamic scheme: %s.", SPH_IMPLEMENTATION);
+ message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours (eta=%f).",
+ kernel_name, p->target_neighbours, p->delta_neighbours,
+ p->eta_neighbours);
+ message("Hydrodynamic integration: CFL parmeter: %.4f.", p->CFL_condition);
+ message(
+ "Hydrodynamic integration: Max change of volume: %.2f (dlog(h)/dt=%f).",
+ powf(expf(p->log_max_h_change), 3.f), p->log_max_h_change);
+}
diff --git a/src/hydro_properties.h b/src/hydro_properties.h
new file mode 100644
index 0000000000..c84252a1dc
--- /dev/null
+++ b/src/hydro_properties.h
@@ -0,0 +1,56 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 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_HYDRO_PROPERTIES
+#define SWIFT_HYDRO_PROPERTIES
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes. */
+#include "const.h"
+#include "parser.h"
+
+/**
+ * @brief Contains all the constants and parameters of the hydro scheme
+ */
+struct hydro_props {
+
+ /* Kernel properties */
+ float eta_neighbours;
+ float target_neighbours;
+ float delta_neighbours;
+
+ /* Kernel properties */
+ int max_smoothing_iterations;
+
+ /* Time integration properties */
+ float CFL_condition;
+ float log_max_h_change;
+
+/* Viscosity parameters */
+#ifdef GADGET_SPH
+ float const_viscosity_alpha;
+#endif
+};
+
+void hydro_props_print(const struct hydro_props *p);
+void hydro_props_init(struct hydro_props *p, const struct swift_params *params);
+
+#endif /* SWIFT_HYDRO_PROPERTIES */
diff --git a/src/kernel_hydro.h b/src/kernel_hydro.h
index 3662510fbc..a9e979151f 100644
--- a/src/kernel_hydro.h
+++ b/src/kernel_hydro.h
@@ -20,6 +20,8 @@
#ifndef SWIFT_KERNEL_HYDRO_H
#define SWIFT_KERNEL_HYDRO_H
+#include <math.h>
+
/* Includes. */
#include "const.h"
#include "error.h"
@@ -143,12 +145,6 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
#define kernel_igamma3 kernel_igamma2 *kernel_igamma
#define kernel_igamma4 kernel_igamma3 *kernel_igamma
-/* Some powers of eta */
-#define kernel_eta3 const_eta_kernel *const_eta_kernel *const_eta_kernel
-
-/* The number of neighbours (i.e. N_ngb) */
-#define kernel_nwneigh 4.0 * M_PI *kernel_gamma3 *kernel_eta3 / 3.0
-
/* Kernel self contribution (i.e. W(0,h)) */
#define kernel_root \
(kernel_coeffs[kernel_degree]) * kernel_constant *kernel_igamma3
diff --git a/src/potentials.h b/src/potentials.h
index 9c687fa4c4..50a8903095 100644
--- a/src/potentials.h
+++ b/src/potentials.h
@@ -31,6 +31,7 @@
#include "const.h"
#include "error.h"
#include "part.h"
+#include "parser.h"
#include "physical_constants.h"
#include "units.h"
diff --git a/src/runner.c b/src/runner.c
index 4e0dad54fa..77153d0608 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -45,6 +45,7 @@
#include "engine.h"
#include "error.h"
#include "gravity.h"
+#include "hydro_properties.h"
#include "hydro.h"
#include "minmax.h"
#include "scheduler.h"
@@ -638,6 +639,13 @@ void runner_doghost(struct runner *r, struct cell *c) {
float h_corr;
const int ti_current = r->e->ti_current;
const double timeBase = r->e->timeBase;
+ const float target_wcount = r->e->hydro_properties->target_neighbours;
+ const float max_wcount =
+ target_wcount + r->e->hydro_properties->delta_neighbours;
+ const float min_wcount =
+ target_wcount - r->e->hydro_properties->delta_neighbours;
+ const int max_smoothing_iter =
+ r->e->hydro_properties->max_smoothing_iterations;
TIMER_TIC;
@@ -654,7 +662,7 @@ void runner_doghost(struct runner *r, struct cell *c) {
for (int k = 0; k < count; k++) pid[k] = k;
/* While there are particles that need to be updated... */
- for (int num_reruns = 0; count > 0 && num_reruns < const_smoothing_max_iter;
+ for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
num_reruns++) {
/* Reset the redo-count. */
@@ -678,7 +686,7 @@ void runner_doghost(struct runner *r, struct cell *c) {
/* Otherwise, compute the smoothing length update (Newton step). */
else {
- h_corr = (kernel_nwneigh - p->density.wcount) / p->density.wcount_dh;
+ h_corr = (target_wcount - p->density.wcount) / p->density.wcount_dh;
/* Truncate to the range [ -p->h/2 , p->h ]. */
h_corr = fminf(h_corr, p->h);
@@ -686,8 +694,7 @@ void runner_doghost(struct runner *r, struct cell *c) {
}
/* Did we get the right number density? */
- if (p->density.wcount > kernel_nwneigh + const_delta_nwneigh ||
- p->density.wcount < kernel_nwneigh - const_delta_nwneigh) {
+ if (p->density.wcount > max_wcount || p->density.wcount < min_wcount) {
/* Ok, correct then */
p->h += h_corr;
@@ -918,7 +925,9 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
struct xpart *const xparts = c->xparts;
struct gpart *const gparts = c->gparts;
const struct external_potential *potential = r->e->external_potential;
+ const struct hydro_props *hydro_properties = r->e->hydro_properties;
const struct phys_const *constants = r->e->physical_constants;
+ const float ln_max_h_change = hydro_properties->log_max_h_change;
const int is_fixdt =
(r->e->policy & engine_policy_fixdt) == engine_policy_fixdt;
@@ -1048,7 +1057,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
} else {
/* Compute the next timestep (hydro condition) */
- const float new_dt_hydro = hydro_compute_timestep(p, xp);
+ const float new_dt_hydro =
+ hydro_compute_timestep(p, xp, hydro_properties);
/* Compute the next timestep (gravity condition) */
float new_dt_grav = FLT_MAX;
@@ -1067,7 +1077,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* Limit change in h */
const float dt_h_change =
- (p->h_dt != 0.0f) ? fabsf(const_ln_max_h_change * p->h / p->h_dt)
+ (p->h_dt != 0.0f) ? fabsf(ln_max_h_change * p->h / p->h_dt)
: FLT_MAX;
new_dt = fminf(new_dt, dt_h_change);
diff --git a/src/swift.h b/src/swift.h
index 4b57044734..7e3116c1de 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -33,6 +33,7 @@
#include "error.h"
#include "gravity.h"
#include "hydro.h"
+#include "hydro_properties.h"
#include "lock.h"
#include "map.h"
#include "multipole.h"
--
GitLab