SPH properties are now read from the parameter file

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

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),

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 \

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"

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 */

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)

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);

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;
 }
 
 /**

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));
 }

src/hydro_properties.c

+/*******************************************************************************
+ * 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);
+}

src/hydro_properties.h

+/*******************************************************************************
+ * 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 */

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

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"
 

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);

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"