Commit cbd001ff authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'master' into timestep_order

parents a84b9c7f c707ab6b
...@@ -368,6 +368,10 @@ int main(int argc, char *argv[]) { ...@@ -368,6 +368,10 @@ int main(int argc, char *argv[]) {
struct hydro_props hydro_properties; struct hydro_props hydro_properties;
if (with_hydro) hydro_props_init(&hydro_properties, params); 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 */ /* Read particles and space information from (GADGET) ICs */
char ICfileName[200] = ""; char ICfileName[200] = "";
parser_get_param_string(params, "InitialConditions:file_name", ICfileName); parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
...@@ -511,7 +515,7 @@ int main(int argc, char *argv[]) { ...@@ -511,7 +515,7 @@ int main(int argc, char *argv[]) {
struct engine e; struct engine e;
engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff, engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff,
engine_policies, talking, reparttype, &us, &prog_const, engine_policies, talking, reparttype, &us, &prog_const,
&hydro_properties, &potential, &cooling_func, &sourceterms); &hydro_properties, &gravity_properties, &potential, &cooling_func, &sourceterms);
if (myrank == 0) { if (myrank == 0) {
clocks_gettime(&toc); clocks_gettime(&toc);
message("engine_init took %.3f %s.", clocks_diff(&tic, &toc), message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
......
...@@ -48,6 +48,14 @@ SPH: ...@@ -48,6 +48,14 @@ SPH:
max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length. 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 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 # Parameters related to the initial conditions
InitialConditions: InitialConditions:
file_name: SedovBlast/sedov.hdf5 # The file to read file_name: SedovBlast/sedov.hdf5 # The file to read
......
...@@ -45,7 +45,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \ ...@@ -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 \ 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 \ 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 \ 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 # Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \ 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 \ ...@@ -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 \ physical_constants.c potential.c hydro_properties.c \
runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \ runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \
statistics.c runner_doiact_vec.c profiler.c dump.c logger.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. # 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 \ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
......
...@@ -3469,6 +3469,7 @@ void engine_init(struct engine *e, struct space *s, ...@@ -3469,6 +3469,7 @@ void engine_init(struct engine *e, struct space *s,
const struct unit_system *internal_units, const struct unit_system *internal_units,
const struct phys_const *physical_constants, const struct phys_const *physical_constants,
const struct hydro_props *hydro, const struct hydro_props *hydro,
const struct gravity_props *gravity,
const struct external_potential *potential, const struct external_potential *potential,
const struct cooling_function_data *cooling_func, const struct cooling_function_data *cooling_func,
struct sourceterms *sourceterms) { struct sourceterms *sourceterms) {
...@@ -3523,6 +3524,7 @@ void engine_init(struct engine *e, struct space *s, ...@@ -3523,6 +3524,7 @@ void engine_init(struct engine *e, struct space *s,
e->wallclock_time = 0.f; e->wallclock_time = 0.f;
e->physical_constants = physical_constants; e->physical_constants = physical_constants;
e->hydro_properties = hydro; e->hydro_properties = hydro;
e->gravity_properties = gravity;
e->external_potential = potential; e->external_potential = potential;
e->cooling_func = cooling_func; e->cooling_func = cooling_func;
e->sourceterms = sourceterms; e->sourceterms = sourceterms;
...@@ -3719,6 +3721,10 @@ void engine_init(struct engine *e, struct space *s, ...@@ -3719,6 +3721,10 @@ void engine_init(struct engine *e, struct space *s,
if (e->policy & engine_policy_hydro) if (e->policy & engine_policy_hydro)
if (e->nodeID == 0) hydro_props_print(e->hydro_properties); 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 */ /* Check we have sensible time bounds */
if (e->timeBegin >= e->timeEnd) if (e->timeBegin >= e->timeEnd)
error( error(
......
...@@ -39,6 +39,7 @@ ...@@ -39,6 +39,7 @@
/* Includes. */ /* Includes. */
#include "clocks.h" #include "clocks.h"
#include "cooling_struct.h" #include "cooling_struct.h"
#include "gravity_properties.h"
#include "parser.h" #include "parser.h"
#include "partition.h" #include "partition.h"
#include "potential.h" #include "potential.h"
...@@ -209,6 +210,9 @@ struct engine { ...@@ -209,6 +210,9 @@ struct engine {
/* Properties of the hydro scheme */ /* Properties of the hydro scheme */
const struct hydro_props *hydro_properties; const struct hydro_props *hydro_properties;
/* Properties of the self-gravity scheme */
const struct gravity_props *gravity_properties;
/* Properties of external gravitational potential */ /* Properties of external gravitational potential */
const struct external_potential *external_potential; const struct external_potential *external_potential;
...@@ -235,6 +239,7 @@ void engine_init(struct engine *e, struct space *s, ...@@ -235,6 +239,7 @@ void engine_init(struct engine *e, struct space *s,
const struct unit_system *internal_units, const struct unit_system *internal_units,
const struct phys_const *physical_constants, const struct phys_const *physical_constants,
const struct hydro_props *hydro, const struct hydro_props *hydro,
const struct gravity_props *gravity,
const struct external_potential *potential, const struct external_potential *potential,
const struct cooling_function_data *cooling, const struct cooling_function_data *cooling,
struct sourceterms *sourceterms); struct sourceterms *sourceterms);
......
/*******************************************************************************
* 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
/*******************************************************************************
* 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 */
...@@ -28,7 +28,6 @@ ...@@ -28,7 +28,6 @@
#endif #endif
/* Local includes. */ /* Local includes. */
#include "const.h"
#include "parser.h" #include "parser.h"
/** /**
...@@ -47,11 +46,6 @@ struct hydro_props { ...@@ -47,11 +46,6 @@ struct hydro_props {
/* Time integration properties */ /* Time integration properties */
float CFL_condition; float CFL_condition;
float log_max_h_change; 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_print(const struct hydro_props *p);
......
...@@ -41,6 +41,7 @@ ...@@ -41,6 +41,7 @@
#include "engine.h" #include "engine.h"
#include "error.h" #include "error.h"
#include "gravity_io.h" #include "gravity_io.h"
#include "gravity_properties.h"
#include "hydro_io.h" #include "hydro_io.h"
#include "hydro_properties.h" #include "hydro_properties.h"
#include "io_properties.h" #include "io_properties.h"
...@@ -750,12 +751,23 @@ void write_output_parallel(struct engine* e, const char* baseName, ...@@ -750,12 +751,23 @@ void write_output_parallel(struct engine* e, const char* baseName,
io_write_code_description(h_file); io_write_code_description(h_file);
/* Print the SPH parameters */ /* Print the SPH parameters */
h_grp = if (e->policy & engine_policy_hydro) {
H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
if (h_grp < 0) error("Error while creating SPH group"); H5P_DEFAULT);
hydro_props_print_snapshot(h_grp, e->hydro_properties); if (h_grp < 0) error("Error while creating SPH group");
writeSPHflavour(h_grp); hydro_props_print_snapshot(h_grp, e->hydro_properties);
H5Gclose(h_grp); 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 */ /* Print the runtime parameters */
h_grp = h_grp =
......
...@@ -41,6 +41,7 @@ ...@@ -41,6 +41,7 @@
#include "engine.h" #include "engine.h"
#include "error.h" #include "error.h"
#include "gravity_io.h" #include "gravity_io.h"
#include "gravity_properties.h"
#include "hydro_io.h" #include "hydro_io.h"
#include "hydro_properties.h" #include "hydro_properties.h"
#include "io_properties.h" #include "io_properties.h"
...@@ -822,12 +823,23 @@ void write_output_serial(struct engine* e, const char* baseName, ...@@ -822,12 +823,23 @@ void write_output_serial(struct engine* e, const char* baseName,
io_write_code_description(h_file); io_write_code_description(h_file);
/* Print the SPH parameters */ /* Print the SPH parameters */
h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT, if (e->policy & engine_policy_hydro) {
H5P_DEFAULT); h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
if (h_grp < 0) error("Error while creating SPH group"); H5P_DEFAULT);
hydro_props_print_snapshot(h_grp, e->hydro_properties); if (h_grp < 0) error("Error while creating SPH group");
writeSPHflavour(h_grp); hydro_props_print_snapshot(h_grp, e->hydro_properties);
H5Gclose(h_grp); 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 */ /* Print the runtime parameters */
h_grp = h_grp =
......
...@@ -40,6 +40,7 @@ ...@@ -40,6 +40,7 @@
#include "engine.h" #include "engine.h"
#include "error.h" #include "error.h"
#include "gravity_io.h" #include "gravity_io.h"
#include "gravity_properties.h"
#include "hydro_io.h" #include "hydro_io.h"
#include "hydro_properties.h" #include "hydro_properties.h"
#include "io_properties.h" #include "io_properties.h"
...@@ -672,12 +673,23 @@ void write_output_single(struct engine* e, const char* baseName, ...@@ -672,12 +673,23 @@ void write_output_single(struct engine* e, const char* baseName,
io_write_code_description(h_file); io_write_code_description(h_file);
/* Print the SPH parameters */ /* Print the SPH parameters */
h_grp = if (e->policy & engine_policy_hydro) {
H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
if (h_grp < 0) error("Error while creating SPH group"); H5P_DEFAULT);
hydro_props_print_snapshot(h_grp, e->hydro_properties); if (h_grp < 0) error("Error while creating SPH group");
writeSPHflavour(h_grp); hydro_props_print_snapshot(h_grp, e->hydro_properties);
H5Gclose(h_grp); 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 */ /* Print the runtime parameters */
h_grp = h_grp =
......
...@@ -35,6 +35,7 @@ ...@@ -35,6 +35,7 @@
#include "engine.h" #include "engine.h"
#include "error.h" #include "error.h"
#include "gravity.h" #include "gravity.h"
#include "gravity_properties.h"
#include "hydro.h" #include "hydro.h"
#include "hydro_properties.h" #include "hydro_properties.h"
#include "lock.h" #include "lock.h"
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment