diff --git a/examples/main.c b/examples/main.c index f7a54363a579eda8d85c50fb5c3b3a9a3bba7df6..f4e6520c800b25f677a992568c7cb113e36b1043 100644 --- a/examples/main.c +++ b/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 diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index 31b948ea4c59de42bce7c3b80b8ae2526b10e0cf..bdcd880c83eb02fdbed8d9d4da8642f3a90fc97b 100644 --- a/examples/parameter_example.yml +++ b/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 diff --git a/src/Makefile.am b/src/Makefile.am index b0615db01c467f243ba3eb0b961370c31a26fe22..3e759d8fd15cfad6b38fe4ed469fb3b745aff786 100644 --- a/src/Makefile.am +++ b/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 \ diff --git a/src/engine.c b/src/engine.c index 6b2a451c9e26a727b67b09eabd38abf28648bb2a..230cfd2163ab00f6bfae8c5ae92134ca38d3663b 100644 --- a/src/engine.c +++ b/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( diff --git a/src/engine.h b/src/engine.h index 1f6a279faf5e147249121b52c3c2ca765d1fed1e..8a719a8952c157e2c567addce2f6fc6491730eca 100644 --- a/src/engine.h +++ b/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); diff --git a/src/gravity_properties.c b/src/gravity_properties.c new file mode 100644 index 0000000000000000000000000000000000000000..6c1ec1f6d9e80d46278892cc7a2081b2cadad923 --- /dev/null +++ b/src/gravity_properties.c @@ -0,0 +1,78 @@ +/******************************************************************************* + * 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 diff --git a/src/gravity_properties.h b/src/gravity_properties.h new file mode 100644 index 0000000000000000000000000000000000000000..6fde91b1d8f2d32a27ec278d5e42f91e5fe930a6 --- /dev/null +++ b/src/gravity_properties.h @@ -0,0 +1,57 @@ +/******************************************************************************* + * 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 */ diff --git a/src/hydro_properties.h b/src/hydro_properties.h index 6b151e2d038fc1ac30e77ad70bc9ef714cec2899..859ae10f18a71ba920f44b8c59e63abd57cd7c92 100644 --- a/src/hydro_properties.h +++ b/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); diff --git a/src/parallel_io.c b/src/parallel_io.c index 1d78410fb9163a4c40e6706a91ce6852f7c35a1d..99e2a50029a58f83c5dd154f19de7265d7924c61 100644 --- a/src/parallel_io.c +++ b/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 = diff --git a/src/serial_io.c b/src/serial_io.c index 28c329070d4995e24a2c375b1967f4e155ea4764..3c2619d7116b8f7ac21913df8e30525b3ab9cb0b 100644 --- a/src/serial_io.c +++ b/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 = diff --git a/src/single_io.c b/src/single_io.c index 264e138728458f6a1d19b57e45f4f8f2e6acdc8f..567b1d941b15dec2d667a2d5eceb169d76f29c45 100644 --- a/src/single_io.c +++ b/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 = diff --git a/src/swift.h b/src/swift.h index c08a4f3209d9eea0fe02ad9112179a0ed7ccae1e..7f1b19b6066c2d55df1cb9101172ae94c9085583 100644 --- a/src/swift.h +++ b/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"