Commit 3b2b2be6 authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'smoothing_length_derivative' into 'master'

Improved multi-timestep SPH

This brings a range of improvement, mostly minor and of code clarity nature to the code. 
It solves:

- Issue #84. Particles are now <128 Bytes in size.
- Issue #85. Smoothing-lengths now evolve more accurately and we implement a time-step criterion more conservative than default Gadget-2 but that is part of EAGLE.
- Issue #88. We now have a 'minimal' hydro implementation that can be used both as a template and an example. The doxygen documentation is based on that version. The documentation is now very extensive. This model is a poor SPH implementation and should only be used to understand the structure of the code and not be seriously run.
- Issue #95. The engine_collect_kick() function now correctly skips empty cells.
- Issue #96. Conserved quantities are now written to a file. We need gravity before this really get useful.
- The description of the SPH flavour in the HDF5 files is now moved to the hydro_io.h files for consistency.


Most of it are straight-forward changes. Note that this does not fix the issue we have in MPI mode.

See merge request !90
parents 2b8930c4 5c966fdf
...@@ -759,7 +759,7 @@ WARN_LOGFILE = ...@@ -759,7 +759,7 @@ WARN_LOGFILE =
# spaces. # spaces.
# Note: If this tag is empty the current directory is searched. # Note: If this tag is empty the current directory is searched.
INPUT = @top_srcdir@ @top_srcdir@/src @top_srcdir@/src/hydro/Default @top_srcdir@/src/gravity/Default INPUT = @top_srcdir@ @top_srcdir@/src @top_srcdir@/src/hydro/Minimal @top_srcdir@/src/gravity/Default
# This tag can be used to specify the character encoding of the source files # This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
......
...@@ -44,11 +44,13 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \ ...@@ -44,11 +44,13 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
kernel.c tools.c part.c kernel.c tools.c part.c
# Include files for distribution, not installation. # Include files for distribution, not installation.
noinst_HEADERS = atomic.h cycle.h error.h inline.h kernel.h vector.h \ noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel.h vector.h \
runner_doiact.h runner_doiact_grav.h units.h intrinsics.h \ runner_doiact.h runner_doiact_grav.h units.h intrinsics.h \
hydro.h hydro_io.h gravity.h \ hydro.h hydro_io.h gravity.h \
gravity/Default/gravity.h gravity/Default/runner_iact_grav.h \ gravity/Default/gravity.h gravity/Default/runner_iact_grav.h \
gravity/Default/gravity_part.h gravity/Default/gravity_part.h
hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h
hydro/Default/hydro.h hydro/Default/hydro_iact.h hydro/Default/hydro_io.h \ hydro/Default/hydro.h hydro/Default/hydro_iact.h hydro/Default/hydro_io.h \
hydro/Default/hydro_debug.h hydro/Default/hydro_part.h \ hydro/Default/hydro_debug.h hydro/Default/hydro_part.h \
hydro/Gadget2/hydro.h hydro/Gadget2/hydro_iact.h hydro/Gadget2/hydro_io.h \ hydro/Gadget2/hydro.h hydro/Gadget2/hydro_iact.h hydro/Gadget2/hydro_io.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_APPROX_MATH_H
#define SWIFT_APPROX_MATH_H
/**
* @brief Approximate version of expf(x) using a 4th order Taylor expansion
*
* The absolute error is of order 10^-6 for -0.2 < x < 0.2.
*
* @param x The number to take the exponential of.
*/
__attribute__((always_inline)) INLINE static float approx_expf(float x) {
return 1.f + x * (1.f + x * (0.5f + x * (1.f / 6.0f + 1.f / 24.0f * x)));
}
#endif /* SWIFT_APPROX_MATH_H */
...@@ -564,6 +564,7 @@ void cell_init_parts(struct cell *c, void *data) { ...@@ -564,6 +564,7 @@ void cell_init_parts(struct cell *c, void *data) {
xp[i].v_full[0] = p[i].v[0]; xp[i].v_full[0] = p[i].v[0];
xp[i].v_full[1] = p[i].v[1]; xp[i].v_full[1] = p[i].v[1];
xp[i].v_full[2] = p[i].v[2]; xp[i].v_full[2] = p[i].v[2];
hydro_first_init_part(&p[i], &xp[i]);
hydro_init_part(&p[i]); hydro_init_part(&p[i]);
hydro_reset_acceleration(&p[i]); hydro_reset_acceleration(&p[i]);
} }
......
...@@ -136,8 +136,8 @@ struct cell { ...@@ -136,8 +136,8 @@ struct cell {
/* Momentum of particles in cell. */ /* Momentum of particles in cell. */
float mom[3], ang[3]; float mom[3], ang[3];
/* Potential and kinetic energy of particles in this cell. */ /* Mass, potential, internal and kinetic energy of particles in this cell. */
double epot, ekin; double mass, e_pot, e_int, e_kin;
/* Number of particles updated in this cell. */ /* Number of particles updated in this cell. */
int updated; int updated;
......
...@@ -267,59 +267,6 @@ void writeAttribute_s(hid_t grp, char* name, const char* str) { ...@@ -267,59 +267,6 @@ void writeAttribute_s(hid_t grp, char* name, const char* str) {
writeStringAttribute(grp, name, str, strlen(str)); writeStringAttribute(grp, name, str, strlen(str));
} }
/* ------------------------------------------------------------------------------------------------
* This part writes the XMF file descriptor enabling a visualisation through
* ParaView
* ------------------------------------------------------------------------------------------------
*/
/**
* @brief Writes the current model of SPH to the file
* @param h_file The (opened) HDF5 file in which to write
*/
void writeSPHflavour(hid_t h_file) {
hid_t h_grpsph = 0;
h_grpsph = H5Gcreate1(h_file, "/SPH", 0);
if (h_grpsph < 0) error("Error while creating SPH group");
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);
#ifdef LEGACY_GADGET2_SPH
writeAttribute_s(h_grpsph, "Thermal Conductivity Model",
"(No treatment) Legacy Gadget-2 as in Springel (2005)");
writeAttribute_s(h_grpsph, "Viscosity Model",
"Legacy Gadget-2 as in Springel (2005)");
writeAttribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
writeAttribute_f(h_grpsph, "Viscosity beta", 3.f);
#else
writeAttribute_s(h_grpsph, "Thermal Conductivity Model",
"Price (2008) without switch");
writeAttribute_f(h_grpsph, "Thermal Conductivity alpha",
const_conductivity_alpha);
writeAttribute_s(h_grpsph, "Viscosity Model",
"Morris & Monaghan (1997), Rosswog, Davies, Thielemann & "
"Piran (2000) with additional Balsara (1995) switch");
writeAttribute_f(h_grpsph, "Viscosity alpha_min", const_viscosity_alpha_min);
writeAttribute_f(h_grpsph, "Viscosity alpha_max", const_viscosity_alpha_max);
writeAttribute_f(h_grpsph, "Viscosity beta", 2.f);
writeAttribute_f(h_grpsph, "Viscosity decay length", const_viscosity_length);
#endif
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);
writeAttribute_s(h_grpsph, "Kernel", kernel_name);
H5Gclose(h_grpsph);
}
/** /**
* @brief Writes the current Unit System * @brief Writes the current Unit System
* @param h_file The (opened) HDF5 file in which to write * @param h_file The (opened) HDF5 file in which to write
...@@ -372,6 +319,12 @@ void writeCodeDescription(hid_t h_file) { ...@@ -372,6 +319,12 @@ void writeCodeDescription(hid_t h_file) {
H5Gclose(h_grpcode); H5Gclose(h_grpcode);
} }
/* ------------------------------------------------------------------------------------------------
* This part writes the XMF file descriptor enabling a visualisation through
* ParaView
* ------------------------------------------------------------------------------------------------
*/
/** /**
* @brief Prepares the XMF file for the new entry * @brief Prepares the XMF file for the new entry
* *
......
...@@ -38,7 +38,7 @@ ...@@ -38,7 +38,7 @@
1.f /* Value taken from (Price,2008), not used in legacy gadget mode */ 1.f /* Value taken from (Price,2008), not used in legacy gadget mode */
/* Time integration constants. */ /* Time integration constants. */
#define const_cfl 0.2f #define const_cfl 0.1f
#define const_ln_max_h_change \ #define const_ln_max_h_change \
0.231111721f /* Particle can't change volume by more than a factor of \ 0.231111721f /* Particle can't change volume by more than a factor of \
2=1.26^3 over one time step */ 2=1.26^3 over one time step */
...@@ -55,7 +55,7 @@ ...@@ -55,7 +55,7 @@
#define const_theta_max \ #define const_theta_max \
0.57735f /* Opening criteria, which is the ratio of the \ 0.57735f /* Opening criteria, which is the ratio of the \
cell distance over the cell width. */ cell distance over the cell width. */
// #define const_G 6.67384e-8f /* Gravitational constant. */
#define const_G 6.672e-8f /* Gravitational constant. */ #define const_G 6.672e-8f /* Gravitational constant. */
#define const_epsilon 0.0014f /* Gravity blending distance. */ #define const_epsilon 0.0014f /* Gravity blending distance. */
#define const_iepsilon 714.285714286f /* Inverse gravity blending distance. */ #define const_iepsilon 714.285714286f /* Inverse gravity blending distance. */
...@@ -66,7 +66,9 @@ ...@@ -66,7 +66,9 @@
#define const_iepsilon6 (const_iepsilon3* const_iepsilon3) #define const_iepsilon6 (const_iepsilon3* const_iepsilon3)
/* SPH variant to use */ /* SPH variant to use */
#define LEGACY_GADGET2_SPH //#define MINIMAL_SPH
#define GADGET2_SPH
//#define DEFAULT_SPH
/* System of units */ /* System of units */
#define const_unit_length_in_cgs 1 /* 3.08567810e16 /\* 1Mpc *\/ */ #define const_unit_length_in_cgs 1 /* 3.08567810e16 /\* 1Mpc *\/ */
......
...@@ -28,10 +28,14 @@ ...@@ -28,10 +28,14 @@
#include "debug.h" #include "debug.h"
/* Import the right hydro definition */ /* Import the right hydro definition */
#ifdef LEGACY_GADGET2_SPH #if defined(MINIMAL_SPH)
#include "./hydro/Minimal/hydro_debug.h"
#elif defined(GADGET2_SPH)
#include "./hydro/Gadget2/hydro_debug.h" #include "./hydro/Gadget2/hydro_debug.h"
#else #elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_debug.h" #include "./hydro/Default/hydro_debug.h"
#else
#error "Invalid choice of SPH variant"
#endif #endif
/** /**
......
...@@ -1616,41 +1616,51 @@ void engine_barrier(struct engine *e, int tid) { ...@@ -1616,41 +1616,51 @@ void engine_barrier(struct engine *e, int tid) {
void engine_collect_kick(struct cell *c) { void engine_collect_kick(struct cell *c) {
int k, updated = 0; int updated = 0;
float t_end_min = FLT_MAX, t_end_max = 0.0f; float t_end_min = FLT_MAX, t_end_max = 0.0f;
double ekin = 0.0, epot = 0.0; double e_kin = 0.0, e_int = 0.0, e_pot = 0.0;
float mom[3] = {0.0f, 0.0f, 0.0f}, ang[3] = {0.0f, 0.0f, 0.0f}; float mom[3] = {0.0f, 0.0f, 0.0f}, ang[3] = {0.0f, 0.0f, 0.0f};
struct cell *cp; struct cell *cp;
/* If I am a super-cell, return immediately. */ /* Skip super-cells (Their values are already set) */
if (c->kick != NULL || c->count == 0) return; if (c->kick != NULL) return;
/* If this cell is not split, I'm in trouble. */ /* Only do something is the cell is non-empty */
if (!c->split) error("Cell has no super-cell."); if (c->count != 0) {
/* Collect the values from the progeny. */ /* If this cell is not split, I'm in trouble. */
for (k = 0; k < 8; k++) if (!c->split) error("Cell has no super-cell.");
if ((cp = c->progeny[k]) != NULL) {
engine_collect_kick(cp); /* Collect the values from the progeny. */
t_end_min = fminf(t_end_min, cp->t_end_min); for (int k = 0; k < 8; k++)
t_end_max = fmaxf(t_end_max, cp->t_end_max); if ((cp = c->progeny[k]) != NULL) {
updated += cp->updated;
ekin += cp->ekin; /* Recurse */
epot += cp->epot; engine_collect_kick(cp);
mom[0] += cp->mom[0];
mom[1] += cp->mom[1]; /* And update */
mom[2] += cp->mom[2]; t_end_min = fminf(t_end_min, cp->t_end_min);
ang[0] += cp->ang[0]; t_end_max = fmaxf(t_end_max, cp->t_end_max);
ang[1] += cp->ang[1]; updated += cp->updated;
ang[2] += cp->ang[2]; e_kin += cp->e_kin;
} e_int += cp->e_int;
e_pot += cp->e_pot;
mom[0] += cp->mom[0];
mom[1] += cp->mom[1];
mom[2] += cp->mom[2];
ang[0] += cp->ang[0];
ang[1] += cp->ang[1];
ang[2] += cp->ang[2];
}
}
/* Store the collected values in the cell. */ /* Store the collected values in the cell. */
c->t_end_min = t_end_min; c->t_end_min = t_end_min;
c->t_end_max = t_end_max; c->t_end_max = t_end_max;
c->updated = updated; c->updated = updated;
c->ekin = ekin; c->e_kin = e_kin;
c->epot = epot; c->e_int = e_int;
c->e_pot = e_pot;
c->mom[0] = mom[0]; c->mom[0] = mom[0];
c->mom[1] = mom[1]; c->mom[1] = mom[1];
c->mom[2] = mom[2]; c->mom[2] = mom[2];
...@@ -1709,12 +1719,12 @@ void engine_init_particles(struct engine *e) { ...@@ -1709,12 +1719,12 @@ void engine_init_particles(struct engine *e) {
message("Initialising particles"); message("Initialising particles");
engine_prepare(e);
/* Make sure all particles are ready to go */ /* Make sure all particles are ready to go */
/* i.e. clean-up any stupid state in the ICs */ /* i.e. clean-up any stupid state in the ICs */
space_map_cells_pre(s, 1, cell_init_parts, NULL); space_map_cells_pre(s, 1, cell_init_parts, NULL);
engine_prepare(e);
engine_marktasks(e); engine_marktasks(e);
// printParticle(e->s->parts, 1000, e->s->nr_parts); // printParticle(e->s->parts, 1000, e->s->nr_parts);
...@@ -1766,8 +1776,9 @@ void engine_init_particles(struct engine *e) { ...@@ -1766,8 +1776,9 @@ void engine_init_particles(struct engine *e) {
engine_launch(e, e->nr_threads, mask, submask); engine_launch(e, e->nr_threads, mask, submask);
TIMER_TOC(timer_runners); TIMER_TOC(timer_runners);
// message("\n0th ENTROPY CONVERSION\n"); // message("\n0th ENTROPY CONVERSION\n")
/* Apply some conversions (e.g. internal energy -> entropy) */
space_map_cells_pre(s, 1, cell_convert_hydro, NULL); space_map_cells_pre(s, 1, cell_convert_hydro, NULL);
// printParticle(e->s->parts, e->s->xparts,1000, e->s->nr_parts); // printParticle(e->s->parts, e->s->xparts,1000, e->s->nr_parts);
...@@ -1787,7 +1798,7 @@ void engine_step(struct engine *e) { ...@@ -1787,7 +1798,7 @@ void engine_step(struct engine *e) {
int k; int k;
int updates = 0; int updates = 0;
float t_end_min = FLT_MAX, t_end_max = 0.f; float t_end_min = FLT_MAX, t_end_max = 0.f;
double epot = 0.0, ekin = 0.0; double e_pot = 0.0, e_int = 0.0, e_kin = 0.0;
float mom[3] = {0.0, 0.0, 0.0}; float mom[3] = {0.0, 0.0, 0.0};
float ang[3] = {0.0, 0.0, 0.0}; float ang[3] = {0.0, 0.0, 0.0};
struct cell *c; struct cell *c;
...@@ -1799,11 +1810,16 @@ void engine_step(struct engine *e) { ...@@ -1799,11 +1810,16 @@ void engine_step(struct engine *e) {
for (k = 0; k < s->nr_cells; k++) for (k = 0; k < s->nr_cells; k++)
if (s->cells[k].nodeID == e->nodeID) { if (s->cells[k].nodeID == e->nodeID) {
c = &s->cells[k]; c = &s->cells[k];
/* Recurse */
engine_collect_kick(c); engine_collect_kick(c);
/* And aggregate */
t_end_min = fminf(t_end_min, c->t_end_min); t_end_min = fminf(t_end_min, c->t_end_min);
t_end_max = fmaxf(t_end_max, c->t_end_max); t_end_max = fmaxf(t_end_max, c->t_end_max);
ekin += c->ekin; e_kin += c->e_kin;
epot += c->epot; e_int += c->e_int;
e_pot += c->e_pot;
updates += c->updated; updates += c->updated;
mom[0] += c->mom[0]; mom[0] += c->mom[0];
mom[1] += c->mom[1]; mom[1] += c->mom[1];
...@@ -1815,7 +1831,7 @@ void engine_step(struct engine *e) { ...@@ -1815,7 +1831,7 @@ void engine_step(struct engine *e) {
/* Aggregate the data from the different nodes. */ /* Aggregate the data from the different nodes. */
#ifdef WITH_MPI #ifdef WITH_MPI
double in[3], out[3]; double in[4], out[4];
out[0] = t_end_min; out[0] = t_end_min;
if (MPI_Allreduce(out, in, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD) != if (MPI_Allreduce(out, in, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD) !=
MPI_SUCCESS) MPI_SUCCESS)
...@@ -1827,20 +1843,16 @@ void engine_step(struct engine *e) { ...@@ -1827,20 +1843,16 @@ void engine_step(struct engine *e) {
error("Failed to aggregate t_end_max."); error("Failed to aggregate t_end_max.");
t_end_max = in[0]; t_end_max = in[0];
out[0] = updates; out[0] = updates;
out[1] = ekin; out[1] = e_kin;
out[2] = epot; out[2] = e_int;
out[3] = e_pot;
if (MPI_Allreduce(out, in, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) != if (MPI_Allreduce(out, in, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
MPI_SUCCESS) MPI_SUCCESS)
error("Failed to aggregate energies."); error("Failed to aggregate energies.");
updates = in[0]; updates = in[0];
ekin = in[1]; e_kin = in[1];
epot = in[2]; e_int = in[2];
/* int nr_parts; e_pot = in[3];
if ( MPI_Allreduce( &s->nr_parts , &nr_parts , 1 , MPI_INT , MPI_SUM ,
MPI_COMM_WORLD ) != MPI_SUCCESS )
error( "Failed to aggregate particle count." );
if ( e->nodeID == 0 )
message( "nr_parts=%i." , nr_parts ); */
#endif #endif
// message("\nDRIFT\n"); // message("\nDRIFT\n");
...@@ -1916,9 +1928,17 @@ if ( e->nodeID == 0 ) ...@@ -1916,9 +1928,17 @@ if ( e->nodeID == 0 )
TIMER_TOC2(timer_step); TIMER_TOC2(timer_step);
if (e->nodeID == 0) { if (e->nodeID == 0) {
/* Print some information to the screen */
printf("%d %f %f %d %.3f\n", e->step, e->time, e->timeStep, updates, printf("%d %f %f %d %.3f\n", e->step, e->time, e->timeStep, updates,
((double)timers[timer_count - 1]) / CPU_TPS * 1000); ((double)timers[timer_count - 1]) / CPU_TPS * 1000);
fflush(stdout); fflush(stdout);
/* Write some energy statistics */
fprintf(e->file_stats, "%d %f %f %f %f %f %f %f %f %f %f %f\n", e->step,
e->time, e_kin, e_int, e_pot, e_kin + e_int + e_pot, mom[0], mom[1],
mom[2], ang[0], ang[1], ang[2]);
fflush(e->file_stats);
} }
// printParticle(e->s->parts, e->s->xparts,1000, e->s->nr_parts); // printParticle(e->s->parts, e->s->xparts,1000, e->s->nr_parts);
...@@ -2211,6 +2231,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads, ...@@ -2211,6 +2231,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
e->timeStep = 0.; e->timeStep = 0.;
e->dt_min = dt_min; e->dt_min = dt_min;
e->dt_max = dt_max; e->dt_max = dt_max;
e->file_stats = NULL;
engine_rank = nodeID; engine_rank = nodeID;
/* Make the space link back to the engine. */ /* Make the space link back to the engine. */
...@@ -2230,6 +2251,14 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads, ...@@ -2230,6 +2251,14 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
#endif #endif
} }
/* Open some files */
if (e->nodeID == 0) {
e->file_stats = fopen("energy.txt", "w");
fprintf(e->file_stats,
"# Step Time E_kin E_int E_pot E_tot "
"p_x p_y p_z ang_x ang_y ang_z\n");
}
/* Print policy */ /* Print policy */
engine_print_policy(e); engine_print_policy(e);
......
...@@ -29,6 +29,7 @@ ...@@ -29,6 +29,7 @@
/* Some standard headers. */ /* Some standard headers. */
#include <pthread.h> #include <pthread.h>
#include <stdio.h>
/* Includes. */ /* Includes. */
#include "lock.h" #include "lock.h"
...@@ -112,8 +113,8 @@ struct engine { ...@@ -112,8 +113,8 @@ struct engine {
/* Time step */ /* Time step */
float timeStep; float timeStep;
/* The system energies from the previous step. */ /* File for statistics */
double ekin, epot; FILE *file_stats;
/* The current step number. */ /* The current step number. */
int step, nullstep; int step, nullstep;
......
...@@ -22,12 +22,17 @@ ...@@ -22,12 +22,17 @@
#include "./const.h" #include "./const.h"
/* Import the right functions */ /* Import the right functions */
#ifdef LEGACY_GADGET2_SPH #if defined(MINIMAL_SPH)
#include "./hydro/Gadget2/hydro.h" #include "./hydro/Minimal/hydro_iact.h"
#include "./hydro/Minimal/hydro.h"
#elif defined(GADGET2_SPH)
#include "./hydro/Gadget2/hydro_iact.h" #include "./hydro/Gadget2/hydro_iact.h"
#else #include "./hydro/Gadget2/hydro.h"
#include "./hydro/Default/hydro.h" #elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_iact.h" #include "./hydro/Default/hydro_iact.h"
#include "./hydro/Default/hydro.h"
#else
#error "Invalid choice of SPH variant"
#endif #endif
#endif #endif /* SWIFT_HYDRO_H */
...@@ -28,19 +28,27 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep( ...@@ -28,19 +28,27 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
struct part* p, struct xpart* xp) { struct part* p, struct xpart* xp) {
/* CFL condition */ /* CFL condition */
float dt_cfl = const_cfl * p->h / p->force.v_sig; const float dt_cfl = 2.f * const_cfl * kernel_gamma * p->h / p->force.v_sig;
/* Limit change in h */
float dt_h_change = (p->force.h_dt != 0.0f)
? fabsf(const_ln_max_h_change * p->h / p->force.h_dt)
: FLT_MAX;
/* Limit change in u */ /* Limit change in u */
float dt_u_change = (p->force.u_dt != 0.0f) const float dt_u_change =
? fabsf(const_max_u_change * p->u / p->force.u_dt) (p->force.u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->force.u_dt)
: FLT_MAX; : FLT_MAX;
return fminf(dt_cfl, fminf(dt_h_change, dt_u_change)); return fminf(dt_cfl, dt_u_change);
}
/**
* @brief Initialises the particles for the first time
*
* This function is called only once just after the ICs have been
* read in to do some conversions.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
__attribute__((always_inline))
INLINE static void hydro_first_init_part(struct part* p, struct xpart* xp) {
} }
/**