Properties of the self-gravity scheme are now read from the YAML file. They...

Properties of the self-gravity scheme are now read from the YAML file. They are printed to stdout and to snapshots.

examples/main.c

@@ -333,10 +333,10 @@ int main(int argc, char *argv[]) {
   MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
 #endif
 
-/* Prepare the domain decomposition scheme */
+  /* Prepare the domain decomposition scheme */
+  enum repartition_type reparttype = REPART_NONE;
 #ifdef WITH_MPI
   struct partition initial_partition;
-  enum repartition_type reparttype;
   partition_init(&initial_partition, &reparttype, params, nr_nodes);
 
   /* Let's report what we did */
@@ -368,6 +368,10 @@ int main(int argc, char *argv[]) {
   struct hydro_props hydro_properties;
   if (with_hydro) hydro_props_init(&hydro_properties, params);
 
+  /* Initialise the gravity properties */
+  struct gravity_props gravity_properties;
+  if (with_self_gravity) gravity_props_init(&gravity_properties, params);
+
   /* Read particles and space information from (GADGET) ICs */
   char ICfileName[200] = "";
   parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
@@ -510,8 +514,9 @@ 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, with_aff,
-              engine_policies, talking, &us, &prog_const, &hydro_properties,
-              &potential, &cooling_func, &sourceterms);
+              engine_policies, talking, reparttype, &us, &prog_const,
+              &hydro_properties, &gravity_properties, &potential, &cooling_func,
+              &sourceterms);
   if (myrank == 0) {
     clocks_gettime(&toc);
     message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
@@ -585,11 +590,6 @@ int main(int argc, char *argv[]) {
   /* Main simulation loop */
   for (int j = 0; !engine_is_done(&e) && e.step != nsteps; j++) {
 
-/* Repartition the space amongst the nodes? */
-#ifdef WITH_MPI
-    if (j % 100 == 2) e.forcerepart = reparttype;
-#endif
-
     /* Reset timers */
     timers_reset(timers_mask_all);
 
@@ -691,6 +691,7 @@ int main(int argc, char *argv[]) {
 #endif
 
   /* Write final output. */
+  engine_drift_all(&e);
   engine_dump_snapshot(&e);
 
 #ifdef WITH_MPI

examples/parameter_example.yml

@@ -48,6 +48,14 @@ SPH:
   max_ghost_iterations:  30       # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
   max_volume_change:     2.       # (Optional) Maximal allowed change of kernel volume over one time-step
 
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                   0.025    # Constant dimensionless multiplier for time integration. 
+  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).
+
+  
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  SedovBlast/sedov.hdf5 # The file to read

src/Makefile.am

@@ -45,7 +45,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     physical_constants.h physical_constants_cgs.h potential.h version.h \
     hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
     sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
-    dump.h logger.h active.h timeline.h xmf.h
+    dump.h logger.h active.h timeline.h xmf.h gravity_properties.h
 
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
@@ -55,7 +55,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     physical_constants.c potential.c hydro_properties.c \
     runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \
     statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
-    part_type.c xmf.c
+    part_type.c xmf.c gravity_properties.c
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \

src/engine.c

@@ -3451,6 +3451,7 @@ void engine_init(struct engine *e, struct space *s,
                  const struct unit_system *internal_units,
                  const struct phys_const *physical_constants,
                  const struct hydro_props *hydro,
+                 const struct gravity_props *gravity,
                  const struct external_potential *potential,
                  const struct cooling_function_data *cooling_func,
                  struct sourceterms *sourceterms) {
@@ -3503,6 +3504,7 @@ void engine_init(struct engine *e, struct space *s,
   e->wallclock_time = 0.f;
   e->physical_constants = physical_constants;
   e->hydro_properties = hydro;
+  e->gravity_properties = gravity;
   e->external_potential = potential;
   e->cooling_func = cooling_func;
   e->sourceterms = sourceterms;
@@ -3699,6 +3701,10 @@ void engine_init(struct engine *e, struct space *s,
   if (e->policy & engine_policy_hydro)
     if (e->nodeID == 0) hydro_props_print(e->hydro_properties);
 
+  /* Print information about the hydro scheme */
+  if (e->policy & engine_policy_self_gravity)
+    if (e->nodeID == 0) gravity_props_print(e->gravity_properties);
+  
   /* Check we have sensible time bounds */
   if (e->timeBegin >= e->timeEnd)
     error(

src/engine.h

@@ -39,6 +39,7 @@
 /* Includes. */
 #include "clocks.h"
 #include "cooling_struct.h"
+#include "gravity_properties.h"
 #include "parser.h"
 #include "partition.h"
 #include "potential.h"
@@ -197,6 +198,9 @@ struct engine {
   /* Properties of the hydro scheme */
   const struct hydro_props *hydro_properties;
 
+  /* Properties of the self-gravity scheme */
+  const struct gravity_props *gravity_properties;
+
   /* Properties of external gravitational potential */
   const struct external_potential *external_potential;
 
@@ -222,6 +226,7 @@ void engine_init(struct engine *e, struct space *s,
                  const struct unit_system *internal_units,
                  const struct phys_const *physical_constants,
                  const struct hydro_props *hydro,
+                 const struct gravity_props *gravity,
                  const struct external_potential *potential,
                  const struct cooling_function_data *cooling,
                  struct sourceterms *sourceterms);

src/gravity_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 "gravity_properties.h"
+
+/* Standard headers */
+#include <float.h>
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "common_io.h"
+#include "dimension.h"
+#include "error.h"
+#include "gravity.h"
+#include "kernel_gravity.h"
+
+#define gravity_props_default_a_smooth 1.25f
+#define gravity_props_default_r_cut 4.5f
+
+void gravity_props_init(struct gravity_props *p,
+                        const struct swift_params *params) {
+
+  /* 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);
+
+  /* Time integration */
+  p->eta = parser_get_param_float(params, "Gravity:eta");
+
+  /* Softening lengths */
+  p->epsilon = parser_get_param_float(params, "Gravity:epsilon");
+}
+
+void gravity_props_print(const struct gravity_props *p) {
+
+  message("Self-gravity scheme: FMM-MM");
+
+  message("Self-gravity time integration: eta=%.4f", p->eta);
+
+  message("Self-gravity softening: epsilon=%.4f", p->epsilon);
+
+  if (p->a_smooth != gravity_props_default_a_smooth)
+    message("Self-gravity smoothing-scale: a_smooth=%f", p->a_smooth);
+
+  if (p->r_cut != gravity_props_default_r_cut)
+    message("Self-gravity MM cut-off: r_cut=%f", p->r_cut);
+}
+
+#if defined(HAVE_HDF5)
+void gravity_props_print_snapshot(hid_t h_grpgrav,
+                                  const struct gravity_props *p) {
+
+  io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta);
+  io_write_attribute_f(h_grpgrav, "Softening", p->epsilon);
+  io_write_attribute_f(h_grpgrav, "MM a_smooth", p->a_smooth);
+  io_write_attribute_f(h_grpgrav, "MM r_cut", p->r_cut);
+}
+#endif

src/gravity_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_GRAVITY_PROPERTIES
+#define SWIFT_GRAVITY_PROPERTIES
+
+/* Config parameters. */
+#include "../config.h"
+
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local includes. */
+#include "parser.h"
+
+/**
+ * @brief Contains all the constants and parameters of the self-gravity scheme
+ */
+struct gravity_props {
+
+  /* Tree-PM parameters */
+  float a_smooth;
+  float r_cut;
+
+  /* Time integration parameters */
+  float eta;
+
+  /* Softening lengths */
+  float epsilon;
+};
+
+void gravity_props_print(const struct gravity_props *p);
+void gravity_props_init(struct gravity_props *p,
+                        const struct swift_params *params);
+
+#if defined(HAVE_HDF5)
+void gravity_props_print_snapshot(hid_t h_grpsph,
+                                  const struct gravity_props *p);
+#endif
+
+#endif /* SWIFT_GRAVITY_PROPERTIES */

src/hydro_properties.h

@@ -28,7 +28,6 @@
 #endif
 
 /* Local includes. */
-#include "const.h"
 #include "parser.h"
 
 /**
@@ -47,11 +46,6 @@ struct hydro_props {
   /* 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);

src/parallel_io.c

@@ -41,6 +41,7 @@
 #include "engine.h"
 #include "error.h"
 #include "gravity_io.h"
+#include "gravity_properties.h"
 #include "hydro_io.h"
 #include "hydro_properties.h"
 #include "io_properties.h"
@@ -750,12 +751,23 @@ void write_output_parallel(struct engine* e, const char* baseName,
   io_write_code_description(h_file);
 
   /* Print the SPH parameters */
-  h_grp =
+  if(e->policy & engine_policy_hydro){
+    h_grp =
       H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-  if (h_grp < 0) error("Error while creating SPH group");
-  hydro_props_print_snapshot(h_grp, e->hydro_properties);
-  writeSPHflavour(h_grp);
-  H5Gclose(h_grp);
+    if (h_grp < 0) error("Error while creating SPH group");
+    hydro_props_print_snapshot(h_grp, e->hydro_properties);
+    writeSPHflavour(h_grp);
+    H5Gclose(h_grp);
+  }
+    
+  /* Print the gravity parameters */
+  if(e->policy & engine_policy_self_gravity) {
+    h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
+		      H5P_DEFAULT);
+    if (h_grp < 0) error("Error while creating gravity group");
+    gravity_props_print_snapshot(h_grp, e->gravity_properties);
+    H5Gclose(h_grp);
+  }
 
   /* Print the runtime parameters */
   h_grp =

src/serial_io.c

@@ -41,6 +41,7 @@
 #include "engine.h"
 #include "error.h"
 #include "gravity_io.h"
+#include "gravity_properties.h"
 #include "hydro_io.h"
 #include "hydro_properties.h"
 #include "io_properties.h"
@@ -822,12 +823,23 @@ void write_output_serial(struct engine* e, const char* baseName,
     io_write_code_description(h_file);
 
     /* Print the SPH parameters */
-    h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
-                      H5P_DEFAULT);
-    if (h_grp < 0) error("Error while creating SPH group");
-    hydro_props_print_snapshot(h_grp, e->hydro_properties);
-    writeSPHflavour(h_grp);
-    H5Gclose(h_grp);
+    if(e->policy & engine_policy_hydro) {
+      h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
+			H5P_DEFAULT);
+      if (h_grp < 0) error("Error while creating SPH group");
+      hydro_props_print_snapshot(h_grp, e->hydro_properties);
+      writeSPHflavour(h_grp);
+      H5Gclose(h_grp);
+    }
+
+    /* Print the gravity parameters */
+    if(e->policy & engine_policy_self_gravity) {
+      h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
+			H5P_DEFAULT);
+      if (h_grp < 0) error("Error while creating gravity group");
+      gravity_props_print_snapshot(h_grp, e->gravity_properties);
+      H5Gclose(h_grp);
+    }
 
     /* Print the runtime parameters */
     h_grp =

src/single_io.c

@@ -40,6 +40,7 @@
 #include "engine.h"
 #include "error.h"
 #include "gravity_io.h"
+#include "gravity_properties.h"
 #include "hydro_io.h"
 #include "hydro_properties.h"
 #include "io_properties.h"
@@ -672,12 +673,23 @@ void write_output_single(struct engine* e, const char* baseName,
   io_write_code_description(h_file);
 
   /* Print the SPH parameters */
-  h_grp =
+  if(e->policy & engine_policy_hydro) {
+    h_grp =
       H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-  if (h_grp < 0) error("Error while creating SPH group");
-  hydro_props_print_snapshot(h_grp, e->hydro_properties);
-  writeSPHflavour(h_grp);
-  H5Gclose(h_grp);
+    if (h_grp < 0) error("Error while creating SPH group");
+    hydro_props_print_snapshot(h_grp, e->hydro_properties);
+    writeSPHflavour(h_grp);
+    H5Gclose(h_grp);
+  }
+
+  /* Print the gravity parameters */
+  if(e->policy & engine_policy_self_gravity) {
+    h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
+		      H5P_DEFAULT);
+    if (h_grp < 0) error("Error while creating gravity group");
+    gravity_props_print_snapshot(h_grp, e->gravity_properties);
+    H5Gclose(h_grp);
+  }
 
   /* Print the runtime parameters */
   h_grp =

src/swift.h

@@ -35,6 +35,7 @@
 #include "engine.h"
 #include "error.h"
 #include "gravity.h"
+#include "gravity_properties.h"
 #include "hydro.h"
 #include "hydro_properties.h"
 #include "lock.h"