diff --git a/examples/CosmoVolume/cosmoVolume.yml b/examples/CosmoVolume/cosmoVolume.yml index 1bc3b3d3a6522710a2f41745a81cac7bae88c575..b857ce4a1e27cb7484148346dafcef2001e67561 100644 --- a/examples/CosmoVolume/cosmoVolume.yml +++ b/examples/CosmoVolume/cosmoVolume.yml @@ -34,11 +34,12 @@ 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.2348 == 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.6 # 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: diff --git a/examples/ExternalPointMass/externalPointMass.yml b/examples/ExternalPointMass/externalPointMass.yml index f834cb7015ff8f94a5595ed754f9d473a4e2d4f8..6d05cb5b43c541490432c8f283b662c9d7c8d184 100644 --- a/examples/ExternalPointMass/externalPointMass.yml +++ b/examples/ExternalPointMass/externalPointMass.yml @@ -34,11 +34,12 @@ 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.2348 == 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: 10. # 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: diff --git a/examples/SedovBlast/sedov.yml b/examples/SedovBlast/sedov.yml index feeee40a1433d5a2f57f9a80c00b18a550a4ebee..b5ad8d855535e70133464fa6a117db37a2c9fb8b 100644 --- a/examples/SedovBlast/sedov.yml +++ b/examples/SedovBlast/sedov.yml @@ -34,11 +34,12 @@ 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.2348 == 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: 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: diff --git a/examples/SodShock/sodShock.yml b/examples/SodShock/sodShock.yml index 6a253a29610768c72dbbf1b9cb0ff35955da477e..7bb1895c15c80c6a4e54248d1d2bf85dea278ce1 100644 --- a/examples/SodShock/sodShock.yml +++ b/examples/SodShock/sodShock.yml @@ -34,11 +34,12 @@ 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.2348 == 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.01 # 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: diff --git a/examples/UniformBox/uniformBox.yml b/examples/UniformBox/uniformBox.yml index 5ff07586f3fe035ec5ccf6474edc732e4483bf10..ce989d22ad8a2310ef86303e1b7b6788fdcc427f 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.2348 == 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 f4dd3d93e30e58bbefa6d243d4e76d32b1194709..e5417c628680bacdc8e896b36639aaa57aed1d28 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/examples/parameter_example.yml b/examples/parameter_example.yml index e629ba78f48c0975ebd4405c071fc70ad3623dce..761e2f67c16766bc341b2bcb47f40ebbc4ea30ad 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -34,11 +34,12 @@ 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.2348 == 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: 3. # Maximal smoothing length allowed (in internal units). + 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: diff --git a/src/Makefile.am b/src/Makefile.am index 57ca1e59742b9a1ce482c56909c1a5684445e32e..f50f5574d16227cc925866596ad1766ee4f5a171 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 793a2e10c23cb0dd355a181e2ad658ce19dc0014..29072173b060a936fe6f4f43c091581d522985b3 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 dd3e34016c1953efffca67e9c7ee1b42f455dec1..6432ef6a9e8107d94c93a32967e291d7e1e4d24d 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,32 +37,20 @@ 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 \ - cell distance over the cell width. */ +#define const_theta_max 0.57735f +#define const_G 6.672e-8f /* Gravitational constant. */ +#define const_epsilon 0.0014f /* Gravity blending distance. */ -#define const_G 6.672e-8f /* Gravitational constant. */ -#define const_epsilon 0.0014f /* Gravity blending distance. */ -#define const_iepsilon 714.285714286f /* Inverse gravity blending distance. */ -#define const_iepsilon2 (const_iepsilon* const_iepsilon) -#define const_iepsilon3 (const_iepsilon2* const_iepsilon) -#define const_iepsilon4 (const_iepsilon2* const_iepsilon2) -#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 diff --git a/src/engine.c b/src/engine.c index d94075558bc91ac68734f67c0758383d7b56e9c2..05b0739c5555a9843f9a72d13b39a7e7831386fa 100644 --- a/src/engine.c +++ b/src/engine.c @@ -2408,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. */ @@ -2415,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 */ @@ -2456,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; @@ -2562,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 f4bc3be00cacf128cbaf58ac5735e17f4fe000e4..9ef7d57599d30aad5e8000e64148812493299d23 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/Default/hydro.h b/src/hydro/Default/hydro.h index 301d957b89702b492bb616e32e6edbe160ab9efb..b4ed2daf385691958d5849b0ac2d651ac6f98d34 100644 --- a/src/hydro/Default/hydro.h +++ b/src/hydro/Default/hydro.h @@ -27,10 +27,14 @@ * */ __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) { + + 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; /* Limit change in u */ const float dt_u_change = diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h index 0e9ad46ddc1d4e8c8d3ffdbf3e81262ec49a7092..09a7f4c7caa6a2a9fcf77c8d83b0ffcab2ef8b54 100644 --- a/src/hydro/Default/hydro_io.h +++ b/src/hydro/Default/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", @@ -123,11 +123,11 @@ void writeSPHflavour(hid_t h_grpsph) { writeAttribute_f(h_grpsph, "Viscosity decay length", const_viscosity_length); /* 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)); writeAttribute_f(h_grpsph, "Maximal Delta u change over dt", const_max_u_change); } diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index 19ab4beb373fa013fc62c8273b2feed35aa36546..9fa1aae9dbb13ac3aa7e2271f30ef57e91ae85ae 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_debug.h b/src/hydro/Gadget2/hydro_debug.h index a4d1f7dd4397ebfc850b582e1ca81fc0d4edb76a..4f00d8e020f9eaac30fe7f272ca1248d7e4eec58 100644 --- a/src/hydro/Gadget2/hydro_debug.h +++ b/src/hydro/Gadget2/hydro_debug.h @@ -23,14 +23,15 @@ __attribute__((always_inline)) "x=[%.3e,%.3e,%.3e], " "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n " "h=%.3e, " - "wcount=%d, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, " + "wcount=%d, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, " + "S=%.3e, " "dS/dt=%.3e, c=%.3e\n" "divV=%.3e, curlV=%.3e, rotV=[%.3e,%.3e,%.3e] \n " "v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n", p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0], xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2], - p->h, (int)p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh, p->rho, - p->force.pressure, p->entropy, p->entropy_dt, p->force.soundspeed, + p->h, (int)p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh, + p->rho, p->force.pressure, p->entropy, p->entropy_dt, p->force.soundspeed, p->div_v, p->force.curl_v, p->density.rot_v[0], p->density.rot_v[1], p->density.rot_v[2], p->force.v_sig, p->h_dt, p->ti_begin, p->ti_end); } diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h index 706fca3258a24eb646a405553249c1c6cf78a8b7..41f3348db10172c1d0d5813debb50365252a413b 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/Minimal/hydro.h b/src/hydro/Minimal/hydro.h index 661f0eee18b46256fdbb693e7b4d76a1375dd92b..01957065cf5b1d89bb839449ac60b9ec965c2c47 100644 --- a/src/hydro/Minimal/hydro.h +++ b/src/hydro/Minimal/hydro.h @@ -27,13 +27,18 @@ * * @param p Pointer to the particle data * @param xp Pointer to the extended particle data + * @param hydro_properties The SPH parameters * */ __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) { + + 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 dt_cfl; } diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h index afe5de83f423e43b4d2480cca1ac3e84d6c549de..43926d526b47fa8f1b92bb862edabefb7417de3c 100644 --- a/src/hydro/Minimal/hydro_io.h +++ b/src/hydro/Minimal/hydro_io.h @@ -104,9 +104,9 @@ 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, "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 */ @@ -115,11 +115,11 @@ void writeSPHflavour(hid_t h_grpsph) { writeAttribute_s(h_grpsph, "Viscosity Model", "No model"); /* 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)); writeAttribute_f(h_grpsph, "Maximal Delta u change over dt", const_max_u_change); } diff --git a/src/hydro_properties.c b/src/hydro_properties.c new file mode 100644 index 0000000000000000000000000000000000000000..3c270b2702edb503c60cdb3cf1029a014d437124 --- /dev/null +++ b/src/hydro_properties.c @@ -0,0 +1,63 @@ +/******************************************************************************* + * 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 <float.h> +#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 = logf(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" + "(max |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 0000000000000000000000000000000000000000..c84252a1dc12f0e5591a7e512fdf4e246f4ab048 --- /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_gravity.h b/src/kernel_gravity.h index 7fd4b061a7e94be01a11b06ad23d9113f579ebb8..cca449efa6dde73eac428d1e20fa5c28156360fa 100644 --- a/src/kernel_gravity.h +++ b/src/kernel_gravity.h @@ -25,9 +25,12 @@ #include "inline.h" #include "vector.h" -/* Gravity kernel stuff - * ----------------------------------------------------------------------------------------------- - */ +#define const_iepsilon (1. / const_epsilon) +#define const_iepsilon2 (const_iepsilon *const_iepsilon) +#define const_iepsilon3 (const_iepsilon2 *const_iepsilon) +#define const_iepsilon4 (const_iepsilon2 *const_iepsilon2) +#define const_iepsilon5 (const_iepsilon3 *const_iepsilon2) +#define const_iepsilon6 (const_iepsilon3 *const_iepsilon3) /* The gravity kernel is defined as a degree 6 polynomial in the distance r. The resulting value should be post-multiplied with r^-3, resulting diff --git a/src/kernel_hydro.h b/src/kernel_hydro.h index 3662510fbcd3817d338d76bc89b1334e3c0cd395..a9e979151fd466ce4538895f52da24af1e48c81b 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 9c687fa4c4a64bcc0ab8cac1a5e495487fa938f1..50a8903095a8a8218e8d02e3f42af07e40be5b43 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 4e0dad54fa3b0725f378bcc7618024206073438b..77153d0608ef086fe592a261e08dd4c701eb3b13 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 4b570447345465684fc51bcd4ca96d07a4b7f518..7e3116c1de8bc8e6cc2f89d0d5cbe9771ffbf33a 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"