diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index 7f8af2fc2679e8bab96f8dd1796794d8dbf6b5fd..8f61a060b37b0e62189160d0a8c61e713cfd3b8f 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -759,7 +759,7 @@ WARN_LOGFILE =
# spaces.
# 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
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
diff --git a/src/Makefile.am b/src/Makefile.am
index a20416c617527b5dfbc2dc324539561d78388f4e..573e51ddc00394a27f5b34bb29252e02a25a9401 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -44,11 +44,13 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
kernel.c tools.c part.c
# 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 \
hydro.h hydro_io.h gravity.h \
gravity/Default/gravity.h gravity/Default/runner_iact_grav.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_debug.h hydro/Default/hydro_part.h \
hydro/Gadget2/hydro.h hydro/Gadget2/hydro_iact.h hydro/Gadget2/hydro_io.h \
diff --git a/src/approx_math.h b/src/approx_math.h
new file mode 100644
index 0000000000000000000000000000000000000000..ef93ea63c383c74caa3eaff65446962872389a35
--- /dev/null
+++ b/src/approx_math.h
@@ -0,0 +1,33 @@
+/*******************************************************************************
+ * 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 */
diff --git a/src/cell.c b/src/cell.c
index f7fe61ba30c7d04362003c495702c9f23a6efbd3..6d84dafdaa4038f6cf30b9198381fab67a8b6b27 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -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[1] = p[i].v[1];
xp[i].v_full[2] = p[i].v[2];
+ hydro_first_init_part(&p[i], &xp[i]);
hydro_init_part(&p[i]);
hydro_reset_acceleration(&p[i]);
}
diff --git a/src/cell.h b/src/cell.h
index e821d0438f96a7dbc20db6da78d00cf48c4c85cc..0ba63d1f037167e49130af64acc05ea1c879b95d 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -136,8 +136,8 @@ struct cell {
/* Momentum of particles in cell. */
float mom[3], ang[3];
- /* Potential and kinetic energy of particles in this cell. */
- double epot, ekin;
+ /* Mass, potential, internal and kinetic energy of particles in this cell. */
+ double mass, e_pot, e_int, e_kin;
/* Number of particles updated in this cell. */
int updated;
diff --git a/src/common_io.c b/src/common_io.c
index b7a9a887e450b3fcc664650688a19910efcfab3d..7059f1e813d6cdaf6afd1e637a826292674a51fe 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -267,59 +267,6 @@ void writeAttribute_s(hid_t grp, char* name, const char* 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
* @param h_file The (opened) HDF5 file in which to write
@@ -372,6 +319,12 @@ void writeCodeDescription(hid_t h_file) {
H5Gclose(h_grpcode);
}
+/* ------------------------------------------------------------------------------------------------
+ * This part writes the XMF file descriptor enabling a visualisation through
+ * ParaView
+ * ------------------------------------------------------------------------------------------------
+ */
+
/**
* @brief Prepares the XMF file for the new entry
*
diff --git a/src/const.h b/src/const.h
index 01326e3a22b3f608dfa9c0cf25a1dffe41697edc..3bd9edff8227a87d040ec7309998364c946307af 100644
--- a/src/const.h
+++ b/src/const.h
@@ -38,7 +38,7 @@
1.f /* Value taken from (Price,2008), not used in legacy gadget mode */
/* Time integration constants. */
-#define const_cfl 0.2f
+#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 */
@@ -55,7 +55,7 @@
#define const_theta_max \
0.57735f /* Opening criteria, which is the ratio of the \
cell distance over the cell width. */
-// #define const_G 6.67384e-8f /* Gravitational constant. */
+
#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. */
@@ -66,7 +66,9 @@
#define const_iepsilon6 (const_iepsilon3* const_iepsilon3)
/* SPH variant to use */
-#define LEGACY_GADGET2_SPH
+//#define MINIMAL_SPH
+#define GADGET2_SPH
+//#define DEFAULT_SPH
/* System of units */
#define const_unit_length_in_cgs 1 /* 3.08567810e16 /\* 1Mpc *\/ */
diff --git a/src/debug.c b/src/debug.c
index 3efca024d10030c3bf4a9c0f9861ec0dc2f07310..2be981175eabe41a8a7e57aec4261a9b38def6ac 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -28,10 +28,14 @@
#include "debug.h"
/* 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"
-#else
+#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_debug.h"
+#else
+#error "Invalid choice of SPH variant"
#endif
/**
diff --git a/src/engine.c b/src/engine.c
index 4fd9fff3ca6eb9346a47f6b4a2b4bc1c5e39b09f..2cca53610c88459e0a076b246645dd4ad998d666 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1616,41 +1616,51 @@ void engine_barrier(struct engine *e, int tid) {
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;
- 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};
struct cell *cp;
- /* If I am a super-cell, return immediately. */
- if (c->kick != NULL || c->count == 0) return;
-
- /* If this cell is not split, I'm in trouble. */
- if (!c->split) error("Cell has no super-cell.");
-
- /* Collect the values from the progeny. */
- for (k = 0; k < 8; k++)
- if ((cp = c->progeny[k]) != NULL) {
- engine_collect_kick(cp);
- t_end_min = fminf(t_end_min, cp->t_end_min);
- t_end_max = fmaxf(t_end_max, cp->t_end_max);
- updated += cp->updated;
- ekin += cp->ekin;
- epot += cp->epot;
- 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];
- }
+ /* Skip super-cells (Their values are already set) */
+ if (c->kick != NULL) return;
+
+ /* Only do something is the cell is non-empty */
+ if (c->count != 0) {
+
+ /* If this cell is not split, I'm in trouble. */
+ if (!c->split) error("Cell has no super-cell.");
+
+ /* Collect the values from the progeny. */
+ for (int k = 0; k < 8; k++)
+ if ((cp = c->progeny[k]) != NULL) {
+
+ /* Recurse */
+ engine_collect_kick(cp);
+
+ /* And update */
+ t_end_min = fminf(t_end_min, cp->t_end_min);
+ t_end_max = fmaxf(t_end_max, cp->t_end_max);
+ updated += cp->updated;
+ 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. */
c->t_end_min = t_end_min;
c->t_end_max = t_end_max;
c->updated = updated;
- c->ekin = ekin;
- c->epot = epot;
+ c->e_kin = e_kin;
+ c->e_int = e_int;
+ c->e_pot = e_pot;
c->mom[0] = mom[0];
c->mom[1] = mom[1];
c->mom[2] = mom[2];
@@ -1709,12 +1719,12 @@ void engine_init_particles(struct engine *e) {
message("Initialising particles");
+ engine_prepare(e);
+
/* Make sure all particles are ready to go */
/* i.e. clean-up any stupid state in the ICs */
space_map_cells_pre(s, 1, cell_init_parts, NULL);
- engine_prepare(e);
-
engine_marktasks(e);
// printParticle(e->s->parts, 1000, e->s->nr_parts);
@@ -1766,8 +1776,9 @@ void engine_init_particles(struct engine *e) {
engine_launch(e, e->nr_threads, mask, submask);
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);
// printParticle(e->s->parts, e->s->xparts,1000, e->s->nr_parts);
@@ -1787,7 +1798,7 @@ void engine_step(struct engine *e) {
int k;
int updates = 0;
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 ang[3] = {0.0, 0.0, 0.0};
struct cell *c;
@@ -1799,11 +1810,16 @@ void engine_step(struct engine *e) {
for (k = 0; k < s->nr_cells; k++)
if (s->cells[k].nodeID == e->nodeID) {
c = &s->cells[k];
+
+ /* Recurse */
engine_collect_kick(c);
+
+ /* And aggregate */
t_end_min = fminf(t_end_min, c->t_end_min);
t_end_max = fmaxf(t_end_max, c->t_end_max);
- ekin += c->ekin;
- epot += c->epot;
+ e_kin += c->e_kin;
+ e_int += c->e_int;
+ e_pot += c->e_pot;
updates += c->updated;
mom[0] += c->mom[0];
mom[1] += c->mom[1];
@@ -1815,7 +1831,7 @@ void engine_step(struct engine *e) {
/* Aggregate the data from the different nodes. */
#ifdef WITH_MPI
- double in[3], out[3];
+ double in[4], out[4];
out[0] = t_end_min;
if (MPI_Allreduce(out, in, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD) !=
MPI_SUCCESS)
@@ -1827,20 +1843,16 @@ void engine_step(struct engine *e) {
error("Failed to aggregate t_end_max.");
t_end_max = in[0];
out[0] = updates;
- out[1] = ekin;
- out[2] = epot;
+ out[1] = e_kin;
+ out[2] = e_int;
+ out[3] = e_pot;
if (MPI_Allreduce(out, in, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate energies.");
updates = in[0];
- ekin = in[1];
- epot = in[2];
-/* int nr_parts;
-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 ); */
+ e_kin = in[1];
+ e_int = in[2];
+ e_pot = in[3];
#endif
// message("\nDRIFT\n");
@@ -1916,9 +1928,17 @@ if ( e->nodeID == 0 )
TIMER_TOC2(timer_step);
if (e->nodeID == 0) {
+
+ /* Print some information to the screen */
printf("%d %f %f %d %.3f\n", e->step, e->time, e->timeStep, updates,
((double)timers[timer_count - 1]) / CPU_TPS * 1000);
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);
@@ -2211,6 +2231,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
e->timeStep = 0.;
e->dt_min = dt_min;
e->dt_max = dt_max;
+ e->file_stats = NULL;
engine_rank = nodeID;
/* 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,
#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 */
engine_print_policy(e);
diff --git a/src/engine.h b/src/engine.h
index 4540a737d39c37133ea515bfbace78143b4ed1d3..46c89a331f1bdf834cf512708503c98fd78b18b9 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -29,6 +29,7 @@
/* Some standard headers. */
#include <pthread.h>
+#include <stdio.h>
/* Includes. */
#include "lock.h"
@@ -112,8 +113,8 @@ struct engine {
/* Time step */
float timeStep;
- /* The system energies from the previous step. */
- double ekin, epot;
+ /* File for statistics */
+ FILE *file_stats;
/* The current step number. */
int step, nullstep;
diff --git a/src/hydro.h b/src/hydro.h
index c57c60d696105c94af70ad1b24fce197163ce698..4b131ea7bb65302aaa69503130f7429710dc221f 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -22,12 +22,17 @@
#include "./const.h"
/* Import the right functions */
-#ifdef LEGACY_GADGET2_SPH
-#include "./hydro/Gadget2/hydro.h"
+#if defined(MINIMAL_SPH)
+#include "./hydro/Minimal/hydro_iact.h"
+#include "./hydro/Minimal/hydro.h"
+#elif defined(GADGET2_SPH)
#include "./hydro/Gadget2/hydro_iact.h"
-#else
-#include "./hydro/Default/hydro.h"
+#include "./hydro/Gadget2/hydro.h"
+#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_iact.h"
+#include "./hydro/Default/hydro.h"
+#else
+#error "Invalid choice of SPH variant"
#endif
-#endif
+#endif /* SWIFT_HYDRO_H */
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 3215390e83db1e5f134768211c52d2dec32cecad..466bc72757be625ab71c90003269823845733291 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -28,19 +28,27 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
struct part* p, struct xpart* xp) {
/* CFL condition */
- float dt_cfl = const_cfl * 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;
+ const float dt_cfl = 2.f * const_cfl * kernel_gamma * p->h / p->force.v_sig;
/* Limit change in u */
- float dt_u_change = (p->force.u_dt != 0.0f)
- ? fabsf(const_max_u_change * p->u / p->force.u_dt)
- : FLT_MAX;
+ const float dt_u_change =
+ (p->force.u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->force.u_dt)
+ : 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) {
}
/**
@@ -97,9 +105,10 @@ __attribute__((always_inline))
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
+ * @param time The current time
*/
-__attribute__((always_inline))
- INLINE static void hydro_prepare_force(struct part* p, struct xpart* xp) {
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+ struct part* p, struct xpart* xp, float time) {
/* Some smoothing length multiples. */
const float h = p->h;
@@ -152,13 +161,13 @@ __attribute__((always_inline))
INLINE static void hydro_reset_acceleration(struct part* p) {
/* Reset the acceleration. */
- p->a[0] = 0.0f;
- p->a[1] = 0.0f;
- p->a[2] = 0.0f;
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
/* Reset the time derivatives. */
p->force.u_dt = 0.0f;
- p->force.h_dt = 0.0f;
+ p->h_dt = 0.0f;
p->force.v_sig = 0.0f;
}
@@ -196,19 +205,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
* @param p The particle to act upon
*/
__attribute__((always_inline))
- INLINE static void hydro_end_force(struct part* p) {
-
- p->force.h_dt *= p->h * 0.333333333f;
-}
+ INLINE static void hydro_end_force(struct part* p) {}
/**
* @brief Kick the additional variables
*
* @param p The particle to act upon
+ * @param xp The particle extended data to act upon
* @param dt The time-step for this kick
+ * @param half_dt The half time-step for this kick
*/
-__attribute__((always_inline))
- INLINE static void hydro_kick_extra(struct part* p, float dt) {}
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+ struct part* p, struct xpart* xp, float dt, float half_dt) { }
/**
* @brief Converts hydro quantity of a particle
@@ -219,3 +227,14 @@ __attribute__((always_inline))
*/
__attribute__((always_inline))
INLINE static void hydro_convert_quantities(struct part* p) {}
+
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline))
+ INLINE static float hydro_get_internal_energy(struct part* p) {
+
+ return p->u;
+}
diff --git a/src/hydro/Default/hydro_debug.h b/src/hydro/Default/hydro_debug.h
index 36eb75ae60e7a49ae798c232377985dd90cff54a..5b1648e222c36dc142362f26ad5188a2af397192 100644
--- a/src/hydro/Default/hydro_debug.h
+++ b/src/hydro/Default/hydro_debug.h
@@ -25,6 +25,7 @@ __attribute__((always_inline))
"h=%.3e, "
"wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, t_begin=%.3e, t_end=%.3e\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[0], p->a[1], p->a[2], 2. * p->h,
- (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->t_begin, p->t_end);
+ 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->mass, p->rho_dh, p->rho, p->t_begin,
+ p->t_end);
}
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index 95342cd852ff359f23fdec755df8df65b41a2ce0..b5b631501b2f9c398cf1f7e5ee32fd5c962ba86e 100644
--- a/src/hydro/Default/hydro_iact.h
+++ b/src/hydro/Default/hydro_iact.h
@@ -432,8 +432,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Use the force, Luke! */
for (k = 0; k < 3; k++) {
f = dx[k] * w;
- pi->a[k] -= mj * f;
- pj->a[k] += mi * f;
+ pi->a_hydro[k] -= mj * f;
+ pj->a_hydro[k] += mi * f;
}
/* Get the time derivative for u. */
@@ -447,8 +447,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->force.u_dt += mi * tc * (pj->u - pi->u);
/* Get the time derivative for h. */
- pi->force.h_dt -= mj * dvdr / rhoj * wi_dr;
- pj->force.h_dt -= mi * dvdr / rhoi * wj_dr;
+ pi->h_dt -= mj * dvdr / rhoj * wi_dr;
+ pj->h_dt -= mi * dvdr / rhoi * wj_dr;
/* Update the signal velocity. */
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
@@ -656,8 +656,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
for (k = 0; k < VEC_SIZE; k++) {
pi[k]->force.u_dt += piu_dt.f[k];
pj[k]->force.u_dt += pju_dt.f[k];
- pi[k]->force.h_dt -= pih_dt.f[k];
- pj[k]->force.h_dt -= pjh_dt.f[k];
+ pi[k]->h_dt -= pih_dt.f[k];
+ pj[k]->h_dt -= pjh_dt.f[k];
pi[k]->force.v_sig = vi_sig.f[k];
pj[k]->force.v_sig = vj_sig.f[k];
for (j = 0; j < 3; j++) {
@@ -747,7 +747,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Use the force, Luke! */
for (k = 0; k < 3; k++) {
f = dx[k] * w;
- pi->a[k] -= mj * f;
+ pi->a_hydro[k] -= mj * f;
}
/* Get the time derivative for u. */
@@ -758,7 +758,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.u_dt += mj * tc * (pi->u - pj->u);
/* Get the time derivative for h. */
- pi->force.h_dt -= mj * dvdr / rhoj * wi_dr;
+ pi->h_dt -= mj * dvdr / rhoj * wi_dr;
/* Update the signal velocity. */
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
@@ -957,7 +957,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
/* Store the forces back on the particles. */
for (k = 0; k < VEC_SIZE; k++) {
pi[k]->force.u_dt += piu_dt.f[k];
- pi[k]->force.h_dt -= pih_dt.f[k];
+ pi[k]->h_dt -= pih_dt.f[k];
pi[k]->force.v_sig = vi_sig.f[k];
pj[k]->force.v_sig = vj_sig.f[k];
for (j = 0; j < 3; j++) pi[k]->a[j] -= pia[j].f[k];
diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h
index fa25a600fd9c928eb27ce93d39d54f1af324df9e..958bf5a1869718b57678246ff3b1985e54145824 100644
--- a/src/hydro/Default/hydro_io.h
+++ b/src/hydro/Default/hydro_io.h
@@ -45,7 +45,7 @@ __attribute__((always_inline)) INLINE static void hydro_read_particles(
COMPULSORY);
readArray(h_grp, "ParticleIDs", ULONGLONG, N, 1, parts, N_total, offset, id,
COMPULSORY);
- readArray(h_grp, "Acceleration", FLOAT, N, 3, parts, N_total, offset, a,
+ readArray(h_grp, "Acceleration", FLOAT, N, 3, parts, N_total, offset, a_hydro,
OPTIONAL);
readArray(h_grp, "Density", FLOAT, N, 1, parts, N_total, offset, rho,
OPTIONAL);
@@ -84,7 +84,43 @@ __attribute__((always_inline)) INLINE static void hydro_write_particles(
writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N, 1, parts,
N_total, mpi_rank, offset, id, us, UNIT_CONV_NO_UNITS);
writeArray(h_grp, fileName, xmfFile, "Acceleration", FLOAT, N, 3, parts,
- N_total, mpi_rank, offset, a, us, UNIT_CONV_ACCELERATION);
+ N_total, mpi_rank, offset, a_hydro, us, UNIT_CONV_ACCELERATION);
writeArray(h_grp, fileName, xmfFile, "Density", FLOAT, N, 1, parts, N_total,
mpi_rank, offset, rho, us, UNIT_CONV_DENSITY);
}
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grpsph The HDF5 group in which to write
+ */
+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);
+
+ /* Viscosity and thermal conduction */
+ 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);
+
+ /* 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, "Maximal Delta u change over dt",
+ const_max_u_change);
+}
diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h
index 6bfe81822f1519b1239428234a72a48868cf170e..b84fc2700b353e1aadc597b6fa7fe1e30676a26d 100644
--- a/src/hydro/Default/hydro_part.h
+++ b/src/hydro/Default/hydro_part.h
@@ -17,10 +17,7 @@
*
******************************************************************************/
-/* Some standard headers. */
-#include <stdlib.h>
-
-/* Extra particle data not needed during the computation. */
+/* Extra particle data not needed during the SPH loops over neighbours. */
struct xpart {
/* Old position, at last tree rebuild. */
@@ -29,9 +26,6 @@ struct xpart {
/* Velocity at the last full step. */
float v_full[3];
- /* Entropy at the half-step. */
- float u_hdt;
-
/* Old density. */
float omega;
@@ -47,11 +41,14 @@ struct part {
float v[3];
/* Particle acceleration. */
- float a[3];
+ float a_hydro[3];
/* Particle cutoff radius. */
float h;
+ /* Change in smoothing length over time. */
+ float h_dt;
+
/* Particle time of beginning of time-step. */
float t_begin;
@@ -101,9 +98,6 @@ struct part {
/* Change in particle energy over time. */
float u_dt;
- /* Change in smoothing length over time. */
- float h_dt;
-
/* Signal velocity */
float v_sig;
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 7ac07f499e1b1997f07dde6361d9f80c4c42ba1e..5194bccc6c3534c9ac896cf97e4cf503b6718273 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -28,17 +28,32 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
struct part* p, struct xpart* xp) {
/* Acceleration */
- float ac = sqrtf(p->a[0] * p->a[0] + p->a[1] * p->a[1] + p->a[2] * p->a[2]);
+ 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
/* CFL condition */
- const float dt_cfl = 2.f * const_cfl * p->h / p->v_sig;
+ const float dt_cfl = 2.f * const_cfl * kernel_gamma * p->h / p->force.v_sig;
return fminf(dt_cfl, dt_accel);
}
+/**
+ * @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) {
+}
+
/**
* @brief Prepares a particle for the density calculation.
*
@@ -54,10 +69,9 @@ __attribute__((always_inline))
p->rho = 0.f;
p->rho_dh = 0.f;
p->div_v = 0.f;
- p->curl_v = 0.f;
- p->rot_v[0] = 0.f;
- p->rot_v[1] = 0.f;
- p->rot_v[2] = 0.f;
+ p->density.rot_v[0] = 0.f;
+ p->density.rot_v[1] = 0.f;
+ p->density.rot_v[2] = 0.f;
}
/**
@@ -86,29 +100,21 @@ __attribute__((always_inline))
/* Finish the calculation by inserting the missing h-factors */
p->rho *= ih * ih2;
p->rho_dh *= ih4;
- p->density.wcount *= (4.0f / 3.0 * M_PI * kernel_gamma3);
- p->density.wcount_dh *= ih * (4.0f / 3.0 * M_PI * kernel_gamma3);
+ p->density.wcount *= (4.0f / 3.0f * M_PI * kernel_gamma3);
+ p->density.wcount_dh *= ih * (4.0f / 3.0f * M_PI * kernel_gamma3);
- /* Compute the derivative term */
- p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh / p->rho);
+ const float irho = 1.f / p->rho;
- /* Finish calculation of the velocity curl */
- p->rot_v[0] *= ih4;
- p->rot_v[1] *= ih4;
- p->rot_v[2] *= ih4;
+ /* Compute the derivative term */
+ p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho);
- /* And compute the norm of the curl */
- p->curl_v = sqrtf(p->rot_v[0] * p->rot_v[0] + p->rot_v[1] * p->rot_v[1] +
- p->rot_v[2] * p->rot_v[2]) /
- p->rho;
+ /* Finish calculation of the velocity curl components */
+ p->density.rot_v[0] *= ih4 * irho;
+ p->density.rot_v[1] *= ih4 * irho;
+ p->density.rot_v[2] *= ih4 * irho;
/* Finish calculation of the velocity divergence */
- p->div_v *= ih4 / p->rho;
-
- /* Compute the pressure */
- const float dt = time - 0.5f * (p->t_begin + p->t_end);
- p->pressure =
- (p->entropy + p->entropy_dt * dt) * powf(p->rho, const_hydro_gamma);
+ p->div_v *= ih4 * irho;
}
/**
@@ -118,38 +124,23 @@ __attribute__((always_inline))
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
+ * @param time The current time
*/
-__attribute__((always_inline))
- INLINE static void hydro_prepare_force(struct part* p, struct xpart* xp) {
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+ struct part* p, struct xpart* xp, float time) {
- /* Some smoothing length multiples. */
- /* const float h = p->h; */
- /* const float ih = 1.0f / h; */
- /* const float ih2 = ih * ih; */
- /* const float ih4 = ih2 * ih2; */
-
- /* /\* Pre-compute some stuff for the balsara switch. *\/ */
- /* const float normDiv_v = fabs(p->density.div_v / p->rho * ih4); */
- /* const float normCurl_v = sqrtf(p->density.curl_v[0] * p->density.curl_v[0]
- * + */
- /* p->density.curl_v[1] * p->density.curl_v[1] +
- */
- /* p->density.curl_v[2] * p->density.curl_v[2]) /
- */
- /* p->rho * ih4; */
-
- /* /\* Compute this particle's sound speed. *\/ */
- /* const float u = p->u; */
- /* const float fc = p->force.c = */
- /* sqrtf(const_hydro_gamma * (const_hydro_gamma - 1.0f) * u); */
-
- /* /\* Compute the P/Omega/rho2. *\/ */
- /* xp->omega = 1.0f + 0.3333333333f * h * p->rho_dh / p->rho; */
- /* p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (p->rho * xp->omega); */
-
- /* /\* Balsara switch *\/ */
- /* p->force.balsara = */
- /* normDiv_v / (normDiv_v + normCurl_v + 0.0001f * fc * ih); */
+ /* Compute the norm of the curl */
+ p->force.curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
+ p->density.rot_v[1] * p->density.rot_v[1] +
+ p->density.rot_v[2] * p->density.rot_v[2]);
+
+ /* Compute the pressure */
+ const float dt = time - 0.5f * (p->t_begin + p->t_end);
+ p->force.pressure =
+ (p->entropy + p->entropy_dt * dt) * powf(p->rho, const_hydro_gamma);
+
+ /* Compute the sound speed */
+ p->force.soundspeed = sqrtf(const_hydro_gamma * p->force.pressure / p->rho);
}
/**
@@ -164,15 +155,17 @@ __attribute__((always_inline))
INLINE static void hydro_reset_acceleration(struct part* p) {
/* Reset the acceleration. */
- p->a[0] = 0.0f;
- p->a[1] = 0.0f;
- p->a[2] = 0.0f;
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+
+ p->h_dt = 0.0f;
/* Reset the time derivatives. */
p->entropy_dt = 0.0f;
/* Reset maximal signal velocity */
- p->v_sig = 0.0f;
+ p->force.v_sig = 0.0f;
}
/**
@@ -186,14 +179,13 @@ __attribute__((always_inline))
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part* p, struct xpart* xp, float t0, float t1) {
- const float dt = t1 - t0;
-
- p->rho *= expf(-p->div_v * dt);
- p->h *= expf(0.33333333f * p->div_v * dt);
-
+ /* Drift the pressure */
const float dt_entr = t1 - 0.5f * (p->t_begin + p->t_end);
- p->pressure =
+ p->force.pressure =
(p->entropy + p->entropy_dt * dt_entr) * powf(p->rho, const_hydro_gamma);
+
+ /* Compute the new sound speed */
+ p->force.soundspeed = sqrtf(const_hydro_gamma * p->force.pressure / p->rho);
}
/**
@@ -214,9 +206,12 @@ __attribute__((always_inline))
* @brief Kick the additional variables
*
* @param p The particle to act upon
+ * @param xp The particle extended data to act upon
+ * @param dt The time-step for this kick
+ * @param half_dt The half time-step for this kick
*/
-__attribute__((always_inline))
- INLINE static void hydro_kick_extra(struct part* p, float dt) {
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+ struct part* p, struct xpart* xp, float dt, float half_dt) {
/* Do not decrease the entropy (temperature) by more than a factor of 2*/
const float entropy_change = p->entropy_dt * dt;
@@ -243,3 +238,15 @@ __attribute__((always_inline))
p->entropy = (const_hydro_gamma - 1.f) * p->entropy *
powf(p->rho, -(const_hydro_gamma - 1.f));
}
+
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline))
+ INLINE static float hydro_get_internal_energy(struct part* p) {
+
+ return p->entropy * powf(p->rho, const_hydro_gamma - 1.f) *
+ (1.f / (const_hydro_gamma - 1.f));
+}
diff --git a/src/hydro/Gadget2/hydro_debug.h b/src/hydro/Gadget2/hydro_debug.h
index 3d17c4f94191cec37521f7f1fd4aa1716c8f5fc2..5ac83e227c3b59595a045ea9fb0474671ee5df92 100644
--- a/src/hydro/Gadget2/hydro_debug.h
+++ b/src/hydro/Gadget2/hydro_debug.h
@@ -24,12 +24,13 @@ __attribute__((always_inline))
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"h=%.3e, "
"wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, "
- "dS/dt=%.3e,\n"
+ "dS/dt=%.3e, c=%.3e\n"
"divV=%.3e, curlV=%.3e, rotV=[%.3e,%.3e,%.3e] \n "
- "v_sig=%e t_begin=%.3e, t_end=%.3e\n",
+ "v_sig=%e dh/dt=%.3e t_begin=%.3e, t_end=%.3e\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[0], p->a[1], p->a[2], 2. * p->h,
- (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->pressure,
- p->entropy, p->entropy_dt, p->div_v, p->curl_v, p->rot_v[0], p->rot_v[1],
- p->rot_v[2], p->v_sig, p->t_begin, p->t_end);
+ 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->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->t_begin, p->t_end);
}
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 31528b63f5b3dc89ce533a3a7acafabee3c91332..f92531bcce30d6c1461bb1366f88bc282a2747e7 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -37,8 +37,6 @@
*errors and interactions
* missed by the Gadget-2 tree-code neighbours search.
*
- * The code uses internal energy instead of entropy as a thermodynamical
- *variable.
*/
/**
@@ -103,13 +101,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
- pi->rot_v[0] += faci * curlvr[0];
- pi->rot_v[1] += faci * curlvr[1];
- pi->rot_v[2] += faci * curlvr[2];
+ pi->density.rot_v[0] += faci * curlvr[0];
+ pi->density.rot_v[1] += faci * curlvr[1];
+ pi->density.rot_v[2] += faci * curlvr[2];
- pj->rot_v[0] += facj * curlvr[0];
- pj->rot_v[1] += facj * curlvr[1];
- pj->rot_v[2] += facj * curlvr[2];
+ pj->density.rot_v[0] += facj * curlvr[0];
+ pj->density.rot_v[1] += facj * curlvr[1];
+ pj->density.rot_v[2] += facj * curlvr[2];
}
/**
@@ -142,9 +140,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
pi->density.wcount += wi;
pi->density.wcount_dh -= u * wi_dx;
- /* const float ih3 = h_inv * h_inv * h_inv; */
- /* const float ih4 = h_inv * h_inv * h_inv * h_inv; */
-
const float fac = mj * wi_dx * ri;
/* Compute dv dot r */
@@ -154,32 +149,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
pi->div_v -= fac * dvdr;
- /* if(pi->id == 515050 && pj->id == 504849) */
- /* message("Interacting with %lld. r=%e hi=%e u=%e W=%e dW/dx=%e dh_drho1=%e
- * dh_drho2=%e\n fac=%e dvdr=%e pj->v=[%.3e,%.3e,%.3e]", */
- /* pj->id, */
- /* r, */
- /* hi, */
- /* u, */
- /* wi * ih3, */
- /* wi_dx * ih4, */
- /* -mj * (3.f * kernel_igamma * wi) * ih4, */
- /* -mj * u * wi_dx * kernel_igamma * ih4, */
- /* fac * ih4, */
- /* dvdr, */
- /* pj->v[0], */
- /* pj->v[1], */
- /* pj->v[2] */
- /* ); */
-
/* Compute dv cross r */
curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
- pi->rot_v[0] += fac * curlvr[0];
- pi->rot_v[1] += fac * curlvr[1];
- pi->rot_v[2] += fac * curlvr[2];
+ pi->density.rot_v[0] += fac * curlvr[0];
+ pi->density.rot_v[1] += fac * curlvr[1];
+ pi->density.rot_v[2] += fac * curlvr[2];
}
/**
@@ -197,12 +174,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float r_inv = 1.0f / r;
/* Get some values in local variables. */
- // const float mi = pi->mass;
+ const float mi = pi->mass;
const float mj = pj->mass;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
- const float pressurei = pi->pressure;
- const float pressurej = pj->pressure;
+ const float pressurei = pi->force.pressure;
+ const float pressurej = pj->force.pressure;
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
@@ -223,70 +200,61 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute sound speeds */
- const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
- const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
- float v_sig = ci + cj;
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
(pi->v[2] - pj->v[2]) * dx[2];
- /* Artificial viscosity term */
- float visc = 0.f;
- if (dvdr < 0.f) {
- const float mu_ij = fac_mu * dvdr * r_inv;
- v_sig -= 3.f * mu_ij;
- const float rho_ij = 0.5f * (rhoi + rhoj);
- const float balsara_i = fabsf(pi->div_v) / (fabsf(pi->div_v) + pi->curl_v +
- 0.0001 * ci / fac_mu / hi);
- const float balsara_j = fabsf(pj->div_v) / (fabsf(pj->div_v) + pj->curl_v +
- 0.0001 * cj / fac_mu / hj);
- visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij / rho_ij *
- (balsara_i + balsara_j);
- }
+ /* Balsara term */
+ const float balsara_i =
+ fabsf(pi->div_v) /
+ (fabsf(pi->div_v) + pi->force.curl_v + 0.0001 * ci / fac_mu / hi);
+ const float balsara_j =
+ fabsf(pj->div_v) /
+ (fabsf(pj->div_v) + pj->force.curl_v + 0.0001 * cj / fac_mu / hj);
+
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = fminf(dvdr, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Signal velocity */
+ const float v_sig = ci + cj - 3.f * mu_ij;
+
+ /* Now construct the full viscosity term */
+ const float rho_ij = 0.5f * (rhoi + rhoj);
+ const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
+ (balsara_i + balsara_j) / rho_ij;
/* Now, convolve with the kernel */
- const float visc_term = 0.5f * mj * visc * (wi_dr + wj_dr) * r_inv;
- const float sph_term =
- mj * (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
+ const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+ const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
/* Eventually got the acceleration */
const float acc = visc_term + sph_term;
- /* //if(pi->id == 1000 && pj->id == 1100) */
- /* if(pi->id == 515050 && pj->id == 504849) */
- /* message("Interacting with %lld. r=%e hi=%e hj=%e dWi/dx=%e dWj/dx=%3e
- * dvdr=%e visc=%e sph=%e", */
- /* pj->id, */
- /* r, */
- /* 2*hi, */
- /* 2*hj, */
- /* wi_dr, */
- /* wj_dr, */
- /* dvdr, */
- /* visc_term, */
- /* sph_term */
- /* ); */
- /* if(pi->id == 1100 && pj->id == 1000) */
- /* message("oO"); */
-
/* Use the force Luke ! */
- pi->a[0] -= acc * dx[0];
- pi->a[1] -= acc * dx[1];
- pi->a[2] -= acc * dx[2];
+ pi->a_hydro[0] -= mj * acc * dx[0];
+ pi->a_hydro[1] -= mj * acc * dx[1];
+ pi->a_hydro[2] -= mj * acc * dx[2];
+
+ pj->a_hydro[0] += mi * acc * dx[0];
+ pj->a_hydro[1] += mi * acc * dx[1];
+ pj->a_hydro[2] += mi * acc * dx[2];
- pj->a[0] += acc * dx[0];
- pj->a[1] += acc * dx[1];
- pj->a[2] += acc * dx[2];
+ /* Get the time derivative for h. */
+ pi->h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+ pj->h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
/* Update the signal velocity. */
- pi->v_sig = fmaxf(pi->v_sig, v_sig);
- pj->v_sig = fmaxf(pj->v_sig, v_sig);
+ pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
+ pj->force.v_sig = fmaxf(pj->force.v_sig, v_sig);
/* Change in entropy */
- pi->entropy_dt += 0.5f * visc_term * dvdr;
- pj->entropy_dt -= 0.5f * visc_term * dvdr;
+ pi->entropy_dt += 0.5f * mj * visc_term * dvdr;
+ pj->entropy_dt -= 0.5f * mi * visc_term * dvdr;
}
/**
@@ -308,8 +276,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float mj = pj->mass;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
- const float pressurei = pi->pressure;
- const float pressurej = pj->pressure;
+ const float pressurei = pi->force.pressure;
+ const float pressurej = pj->force.pressure;
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
@@ -330,47 +298,54 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute sound speeds */
- const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
- const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
- float v_sig = ci + cj;
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
(pi->v[2] - pj->v[2]) * dx[2];
- /* Artificial viscosity term */
- float visc = 0.f;
- if (dvdr < 0.f) {
- const float mu_ij = fac_mu * dvdr * r_inv;
- v_sig -= 3.f * mu_ij;
- const float rho_ij = 0.5f * (rhoi + rhoj);
- const float balsara_i = fabsf(pi->div_v) / (fabsf(pi->div_v) + pi->curl_v +
- 0.0001 * ci / fac_mu / hi);
- const float balsara_j = fabsf(pj->div_v) / (fabsf(pj->div_v) + pj->curl_v +
- 0.0001 * cj / fac_mu / hj);
- visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij / rho_ij *
- (balsara_i + balsara_j);
- }
+ /* Balsara term */
+ const float balsara_i =
+ fabsf(pi->div_v) /
+ (fabsf(pi->div_v) + pi->force.curl_v + 0.0001 * ci / fac_mu / hi);
+ const float balsara_j =
+ fabsf(pj->div_v) /
+ (fabsf(pj->div_v) + pj->force.curl_v + 0.0001 * cj / fac_mu / hj);
+
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = fminf(dvdr, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Signal velocity */
+ const float v_sig = ci + cj - 3.f * mu_ij;
+
+ /* Now construct the full viscosity term */
+ const float rho_ij = 0.5f * (rhoi + rhoj);
+ const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
+ (balsara_i + balsara_j) / rho_ij;
/* Now, convolve with the kernel */
- const float visc_term = 0.5f * mj * visc * (wi_dr + wj_dr) * r_inv;
- const float sph_term =
- mj * (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
+ const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+ const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
/* Eventually got the acceleration */
const float acc = visc_term + sph_term;
/* Use the force Luke ! */
- pi->a[0] -= acc * dx[0];
- pi->a[1] -= acc * dx[1];
- pi->a[2] -= acc * dx[2];
+ pi->a_hydro[0] -= mj * acc * dx[0];
+ pi->a_hydro[1] -= mj * acc * dx[1];
+ pi->a_hydro[2] -= mj * acc * dx[2];
+
+ /* Get the time derivative for h. */
+ pi->h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
/* Update the signal velocity. */
- pi->v_sig = fmaxf(pi->v_sig, v_sig);
+ pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
/* Change in entropy */
- pi->entropy_dt += 0.5f * visc_term * dvdr;
+ pi->entropy_dt += 0.5f * mj * visc_term * dvdr;
}
#endif /* SWIFT_RUNNER_IACT_LEGACY_H */
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 5e01d67fddcc81d766e8aa59eee2d6368378170c..17c3d3013644c3572f3c26fc3e270b1c1bc465ed 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -45,7 +45,7 @@ __attribute__((always_inline)) INLINE static void hydro_read_particles(
entropy, COMPULSORY);
readArray(h_grp, "ParticleIDs", ULONGLONG, N, 1, parts, N_total, offset, id,
COMPULSORY);
- readArray(h_grp, "Acceleration", FLOAT, N, 3, parts, N_total, offset, a,
+ readArray(h_grp, "Acceleration", FLOAT, N, 3, parts, N_total, offset, a_hydro,
OPTIONAL);
readArray(h_grp, "Density", FLOAT, N, 1, parts, N_total, offset, rho,
OPTIONAL);
@@ -85,7 +85,36 @@ __attribute__((always_inline)) INLINE static void hydro_write_particles(
writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N, 1, parts,
N_total, mpi_rank, offset, id, us, UNIT_CONV_NO_UNITS);
writeArray(h_grp, fileName, xmfFile, "Acceleration", FLOAT, N, 3, parts,
- N_total, mpi_rank, offset, a, us, UNIT_CONV_ACCELERATION);
+ N_total, mpi_rank, offset, a_hydro, us, UNIT_CONV_ACCELERATION);
writeArray(h_grp, fileName, xmfFile, "Density", FLOAT, N, 1, parts, N_total,
mpi_rank, offset, rho, us, UNIT_CONV_DENSITY);
}
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grpsph The HDF5 group in which to write
+ */
+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);
+
+ /* Viscosity and thermal conduction */
+ 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);
+
+ /* 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));
+}
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 8ad6a6ef24add67f6fcfea0059738b44b3f20175..41898513aadddf7a7f7f5f6acb7516af6ecad1f4 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -17,10 +17,7 @@
*
******************************************************************************/
-/* Some standard headers. */
-#include <stdlib.h>
-
-/* Extra particle data not needed during the computation. */
+/* Extra particle data not needed during the SPH loops over neighbours. */
struct xpart {
/* Old position, at last tree rebuild. */
@@ -29,12 +26,6 @@ struct xpart {
/* Velocity at the last full step. */
float v_full[3];
- /* Entropy at the half-step. */
- float u_hdt;
-
- /* Old density. */
- float omega;
-
} __attribute__((aligned(xpart_align)));
/* Data of a single particle. */
@@ -47,30 +38,20 @@ struct part {
float v[3];
/* Particle acceleration. */
- float a[3];
+ float a_hydro[3];
/* Particle cutoff radius. */
float h;
+ /* Time derivative of the smoothing length */
+ float h_dt;
+
/* Particle time of beginning of time-step. */
float t_begin;
/* Particle time of end of time-step. */
float t_end;
- struct {
-
- /* Number of neighbours */
- float wcount;
-
- /* Number of neighbours spatial derivative */
- float wcount_dh;
-
- } density;
-
- /* Particle entropy. */
- float entropy;
-
/* Particle density. */
float rho;
@@ -78,24 +59,49 @@ struct part {
*/
float rho_dh;
+ /* Particle entropy. */
+ float entropy;
+
+ /* Entropy time derivative */
+ float entropy_dt;
+
/* Particle mass. */
float mass;
- /* Particle pressure */
- float pressure;
+ union {
- /* Entropy time derivative */
- float entropy_dt;
+ struct {
- /* Velocity divergence */
- float div_v;
+ /* Number of neighbours */
+ float wcount;
+
+ /* Number of neighbours spatial derivative */
+ float wcount_dh;
+
+ /* Velocity curl components */
+ float rot_v[3];
+
+ } density;
- /* Velocity curl */
- float curl_v;
- float rot_v[3];
+ struct {
- /* Signal velocity */
- float v_sig;
+ /* Velocity curl norm*/
+ float curl_v;
+
+ /* Signal velocity */
+ float v_sig;
+
+ /* Particle pressure */
+ float pressure;
+
+ /* Particle sound speed */
+ float soundspeed;
+
+ } force;
+ };
+
+ /* Velocity divergence */
+ float div_v;
/* Particle ID. */
unsigned long long id;
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
new file mode 100644
index 0000000000000000000000000000000000000000..3585fc3c295be2d449adee5425d61c6024fa7caf
--- /dev/null
+++ b/src/hydro/Minimal/hydro.h
@@ -0,0 +1,235 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+
+#include "approx_math.h"
+
+/**
+ * @brief Computes the hydro time-step of a given particle
+ *
+ * This function returns the time-step of a particle given its hydro-dynamical
+ * state. A typical time-step calculation would be the use of the CFL condition.
+ *
+ * @param p Pointer to the particle data
+ * @param xp Pointer to the extended particle data
+ *
+ */
+__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
+ struct part* p, struct xpart* xp) {
+
+ /* CFL condition */
+ const float dt_cfl = 2.f * const_cfl * kernel_gamma * p->h / p->force.v_sig;
+
+ return dt_cfl;
+}
+
+/**
+ * @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 or assignments between the particle
+ * and extended particle fields.
+ *
+ * @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) {
+
+ xp->u_full = p->u;
+}
+
+/**
+ * @brief Prepares a particle for the density calculation.
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the various density loop over neighbours. Typically, all fields of the
+ * density sub-structure of a particle get zeroed in here.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline))
+ INLINE static void hydro_init_part(struct part* p) {
+ p->density.wcount = 0.f;
+ p->density.wcount_dh = 0.f;
+ p->rho = 0.f;
+ p->rho_dh = 0.f;
+}
+
+/**
+ * @brief Finishes the density calculation.
+ *
+ * Multiplies the density and number of neighbours by the appropiate constants
+ * and add the self-contribution term.
+ * Additional quantities such as velocity gradients will also get the final
+ *terms
+ * added to them here.
+ *
+ * @param p The particle to act upon
+ * @param time The current time
+ */
+__attribute__((always_inline))
+ INLINE static void hydro_end_density(struct part* p, float time) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float ih = 1.0f / h;
+ const float ih2 = ih * ih;
+ const float ih4 = ih2 * ih2;
+
+ /* Final operation on the density (add self-contribution). */
+ p->rho += p->mass * kernel_root;
+ p->rho_dh -= 3.0f * p->mass * kernel_root * kernel_igamma;
+ p->density.wcount += kernel_root;
+
+ /* Finish the calculation by inserting the missing h-factors */
+ p->rho *= ih * ih2;
+ p->rho_dh *= ih4;
+ p->density.wcount *= (4.0f / 3.0f * M_PI * kernel_gamma3);
+ p->density.wcount_dh *= ih * (4.0f / 3.0f * M_PI * kernel_gamma3);
+}
+
+/**
+ * @brief Prepare a particle for the force calculation.
+ *
+ * This function is called in the ghost task to convert some quantities coming
+ * from the density loop over neighbours into quantities ready to be used in the
+ * force loop over neighbours. Quantities are typically read from the density
+ * sub-structure and written to the force sub-structure.
+ * Examples of calculations done here include the calculation of viscosity term
+ * constants, thermal conduction terms, hydro conversions, etc.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param time The current time
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+ struct part* p, struct xpart* xp, float time) {
+
+ p->force.pressure = p->rho * p->u * (const_hydro_gamma - 1.f);
+}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * Resets all hydro acceleration and time derivative fields in preparation
+ * for the sums taking place in the various force tasks.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline))
+ INLINE static void hydro_reset_acceleration(struct part* p) {
+
+ /* Reset the acceleration. */
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+
+ /* Reset the time derivatives. */
+ p->u_dt = 0.0f;
+ p->h_dt = 0.0f;
+ p->force.v_sig = 0.0f;
+}
+
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * Additional hydrodynamic quantites are drifted forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * @param p The particle
+ * @param xp The extended data of the particle
+ * @param t0 The time at the start of the drift
+ * @param t1 The time at the end of the drift
+ */
+__attribute__((always_inline)) INLINE static void hydro_predict_extra(
+ struct part* p, struct xpart* xp, float t0, float t1) {
+
+ const float dt = t1 - t0;
+
+ /* Predict internal energy */
+ const float w = p->u_dt / p->u * dt;
+ if (fabsf(w) < 0.2f)
+ p->u *= approx_expf(w); /* 4th order expansion of exp(w) */
+ else
+ p->u *= expf(w);
+
+ /* Need to recompute the pressure as well */
+ p->force.pressure = p->rho * p->u * (const_hydro_gamma - 1.f);
+}
+
+/**
+ * @brief Finishes the force calculation.
+ *
+ * Multiplies the force and accelerations by the appropiate constants
+ * and add the self-contribution term. In most cases, there is nothing
+ * to do here.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline))
+ INLINE static void hydro_end_force(struct part* p) {}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * Additional hydrodynamic quantites are kicked forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * @param p The particle to act upon
+ * @param xp The particle extended data to act upon
+ * @param dt The time-step for this kick
+ * @param half_dt The half time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+ struct part* p, struct xpart* xp, float dt, float half_dt) {
+
+ /* Kick in momentum space */
+ xp->u_full += p->u_dt * dt;
+
+ /* Get the predicted internal energy */
+ p->u = xp->u_full - half_dt * p->u_dt;
+}
+
+/**
+ * @brief Converts hydro quantity of a particle at the start of a run
+ *
+ * This function is called once at the end of the engine_init_particle()
+ * routine (at the start of a calculation) after the densities of
+ * particles have been computed.
+ * This can be used to convert internal energy into entropy for instance.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline))
+ INLINE static void hydro_convert_quantities(struct part* p) {}
+
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline))
+ INLINE static float hydro_get_internal_energy(struct part* p) {
+
+ return p->u;
+}
diff --git a/src/hydro/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..b208d5b867779aaaa204d8783b75681156506fe7
--- /dev/null
+++ b/src/hydro/Minimal/hydro_debug.h
@@ -0,0 +1,33 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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/>.
+ *
+ ******************************************************************************/
+
+__attribute__((always_inline))
+ INLINE static void hydro_debug_particle(struct part* p, struct xpart* xp) {
+ printf(
+ "x=[%.3e,%.3e,%.3e], "
+ "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
+ "u_full=%.3e, u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
+ "h=%.3e, dh/dt=%.3e "
+ "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, t_begin=%.3e, t_end=%.3e\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],
+ xp->u_full, p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h,
+ p->h_dt, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->t_begin,
+ p->t_end);
+}
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..6afb9d8d38a4fc7f1d38b7286720ddb7f3c51ab4
--- /dev/null
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -0,0 +1,248 @@
+/*******************************************************************************
+ * 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_RUNNER_IACT_H
+#define SWIFT_RUNNER_IACT_H
+
+/* Includes. */
+#include "const.h"
+#include "kernel.h"
+#include "part.h"
+#include "vector.h"
+
+/**
+ * @brief Minimal conservative implementation of SPH
+ *
+ * The thermal variable is the internal energy (u). No viscosity nor
+ * thermal conduction terms are implemented.
+ */
+
+/**
+ * @brief Density loop
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_density(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ float r = sqrtf(r2);
+ float xi, xj;
+ float h_inv;
+ float wi, wj, wi_dx, wj_dx;
+ float mi, mj;
+
+ /* Get the masses. */
+ mi = pi->mass;
+ mj = pj->mass;
+
+ /* Compute density of pi. */
+ h_inv = 1.0 / hi;
+ xi = r * h_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ pi->rho += mj * wi;
+ pi->rho_dh -= mj * (3.0 * wi + xi * wi_dx);
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= xi * wi_dx;
+
+ /* Compute density of pj. */
+ h_inv = 1.f / hj;
+ xj = r * h_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ pj->rho += mi * wj;
+ pj->rho_dh -= mi * (3.0 * wj + xj * wj_dx);
+ pj->density.wcount += wj;
+ pj->density.wcount_dh -= xj * wj_dx;
+}
+
+/**
+ * @brief Density loop (non-symmetric version)
+ */
+
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ float r;
+ float xi;
+ float h_inv;
+ float wi, wi_dx;
+ float mj;
+
+ /* Get the masses. */
+ mj = pj->mass;
+
+ /* Get r and r inverse. */
+ r = sqrtf(r2);
+
+ h_inv = 1.f / hi;
+ xi = r * h_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ pi->rho += mj * wi;
+ pi->rho_dh -= mj * (3.0 * wi + xi * wi_dx);
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= xi * wi_dx;
+}
+
+/**
+ * @brief Force loop
+ */
+
+__attribute__((always_inline)) INLINE static void runner_iact_force(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
+
+ /* Recover some data */
+ const float mi = pi->mass;
+ const float mj = pj->mass;
+ const float rhoi = pi->rho;
+ const float rhoj = pj->rho;
+ const float pressurei = pi->force.pressure;
+ const float pressurej = pj->force.pressure;
+
+ /* Get the kernel for hi. */
+ const float hi_inv = 1.0f / hi;
+ const float hi2_inv = hi_inv * hi_inv;
+ const float xi = r * hi_inv;
+ float wi, wi_dx;
+ kernel_deval(xi, &wi, &wi_dx);
+ const float wi_dr = hi2_inv * hi2_inv * wi_dx;
+
+ /* Get the kernel for hj. */
+ const float hj_inv = 1.0f / hj;
+ const float hj2_inv = hj_inv * hj_inv;
+ const float xj = r * hj_inv;
+ float wj, wj_dx;
+ kernel_deval(xj, &wj, &wj_dx);
+ const float wj_dr = hj2_inv * hj2_inv * wj_dx;
+
+ /* Compute gradient terms */
+ const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
+ const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
+
+ /* Compute dv dot r. */
+ float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
+ dvdr *= r_inv;
+
+ /* Compute the relative velocity. (This is 0 if the particles move away from
+ * each other and negative otherwise) */
+ const float omega_ij = fminf(dvdr, 0.f);
+
+ /* Compute sound speeds */
+ const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
+ const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
+ float v_sig = ci + cj + 3.f * omega_ij;
+
+ /* SPH acceleration term */
+ const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
+
+ /* Use the force Luke ! */
+ pi->a_hydro[0] -= mj * sph_term * dx[0];
+ pi->a_hydro[1] -= mj * sph_term * dx[1];
+ pi->a_hydro[2] -= mj * sph_term * dx[2];
+
+ pj->a_hydro[0] += mi * sph_term * dx[0];
+ pj->a_hydro[1] += mi * sph_term * dx[1];
+ pj->a_hydro[2] += mi * sph_term * dx[2];
+
+ /* Get the time derivative for u. */
+ pi->u_dt += P_over_rho_i * mj * dvdr * wi_dr;
+ pj->u_dt += P_over_rho_j * mi * dvdr * wj_dr;
+
+ /* Get the time derivative for h. */
+ pi->h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+ pj->h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
+
+ /* Update the signal velocity. */
+ pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
+ pj->force.v_sig = fmaxf(pj->force.v_sig, v_sig);
+}
+
+/**
+ * @brief Force loop (non-symmetric version)
+ */
+
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
+
+ /* Recover some data */
+ // const float mi = pi->mass;
+ const float mj = pj->mass;
+ const float rhoi = pi->rho;
+ const float rhoj = pj->rho;
+ const float pressurei = pi->force.pressure;
+ const float pressurej = pj->force.pressure;
+
+ /* Get the kernel for hi. */
+ const float hi_inv = 1.0f / hi;
+ const float hi2_inv = hi_inv * hi_inv;
+ const float xi = r * hi_inv;
+ float wi, wi_dx;
+ kernel_deval(xi, &wi, &wi_dx);
+ const float wi_dr = hi2_inv * hi2_inv * wi_dx;
+
+ /* Get the kernel for hj. */
+ const float hj_inv = 1.0f / hj;
+ const float hj2_inv = hj_inv * hj_inv;
+ const float xj = r * hj_inv;
+ float wj, wj_dx;
+ kernel_deval(xj, &wj, &wj_dx);
+ const float wj_dr = hj2_inv * hj2_inv * wj_dx;
+
+ /* Compute gradient terms */
+ const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
+ const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
+
+ /* Compute dv dot r. */
+ float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
+ dvdr *= r_inv;
+
+ /* Compute the relative velocity. (This is 0 if the particles move away from
+ * each other and negative otherwise) */
+ const float omega_ij = fminf(dvdr, 0.f);
+
+ /* Compute sound speeds */
+ const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
+ const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
+ float v_sig = ci + cj + 3.f * omega_ij;
+
+ /* SPH acceleration term */
+ const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
+
+ /* Use the force Luke ! */
+ pi->a_hydro[0] -= mj * sph_term * dx[0];
+ pi->a_hydro[1] -= mj * sph_term * dx[1];
+ pi->a_hydro[2] -= mj * sph_term * dx[2];
+
+ /* Get the time derivative for u. */
+ pi->u_dt += P_over_rho_i * mj * dvdr * wi_dr;
+
+ /* Get the time derivative for h. */
+ pi->h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+
+ /* Update the signal velocity. */
+ pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
+}
+
+#endif /* SWIFT_RUNNER_IACT_H */
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..2c56fb489ab84ca7c30426b54cf95e26e3821084
--- /dev/null
+++ b/src/hydro/Minimal/hydro_io.h
@@ -0,0 +1,118 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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/>.
+ *
+ ******************************************************************************/
+
+/**
+ * @brief Reads the different particles to the HDF5 file
+ *
+ * @param h_grp The HDF5 group in which to read the arrays.
+ * @param N The number of particles on that MPI rank.
+ * @param N_total The total number of particles (only used in MPI mode)
+ * @param offset The offset of the particles for this MPI rank (only used in MPI
+ *mode)
+ * @param parts The particle array
+ *
+ */
+__attribute__((always_inline)) INLINE static void hydro_read_particles(
+ hid_t h_grp, int N, long long N_total, long long offset,
+ struct part* parts) {
+
+ /* Read arrays */
+ readArray(h_grp, "Coordinates", DOUBLE, N, 3, parts, N_total, offset, x,
+ COMPULSORY);
+ readArray(h_grp, "Velocities", FLOAT, N, 3, parts, N_total, offset, v,
+ COMPULSORY);
+ readArray(h_grp, "Masses", FLOAT, N, 1, parts, N_total, offset, mass,
+ COMPULSORY);
+ readArray(h_grp, "SmoothingLength", FLOAT, N, 1, parts, N_total, offset, h,
+ COMPULSORY);
+ readArray(h_grp, "InternalEnergy", FLOAT, N, 1, parts, N_total, offset, u,
+ COMPULSORY);
+ readArray(h_grp, "ParticleIDs", ULONGLONG, N, 1, parts, N_total, offset, id,
+ COMPULSORY);
+ readArray(h_grp, "Acceleration", FLOAT, N, 3, parts, N_total, offset, a_hydro,
+ OPTIONAL);
+ readArray(h_grp, "Density", FLOAT, N, 1, parts, N_total, offset, rho,
+ OPTIONAL);
+}
+
+/**
+ * @brief Writes the different particles to the HDF5 file
+ *
+ * @param h_grp The HDF5 group in which to write the arrays.
+ * @param fileName The name of the file (unsued in MPI mode).
+ * @param xmfFile The XMF file to write to (unused in MPI mode).
+ * @param N The number of particles on that MPI rank.
+ * @param N_total The total number of particles (only used in MPI mode)
+ * @param mpi_rank The MPI rank of this node (only used in MPI mode)
+ * @param offset The offset of the particles for this MPI rank (only used in MPI
+ *mode)
+ * @param parts The particle array
+ * @param us The unit system to use
+ *
+ */
+__attribute__((always_inline)) INLINE static void hydro_write_particles(
+ hid_t h_grp, char* fileName, FILE* xmfFile, int N, long long N_total,
+ int mpi_rank, long long offset, struct part* parts, struct UnitSystem* us) {
+
+ /* Write arrays */
+ writeArray(h_grp, fileName, xmfFile, "Coordinates", DOUBLE, N, 3, parts,
+ N_total, mpi_rank, offset, x, us, UNIT_CONV_LENGTH);
+ writeArray(h_grp, fileName, xmfFile, "Velocities", FLOAT, N, 3, parts,
+ N_total, mpi_rank, offset, v, us, UNIT_CONV_SPEED);
+ writeArray(h_grp, fileName, xmfFile, "Masses", FLOAT, N, 1, parts, N_total,
+ mpi_rank, offset, mass, us, UNIT_CONV_MASS);
+ writeArray(h_grp, fileName, xmfFile, "SmoothingLength", FLOAT, N, 1, parts,
+ N_total, mpi_rank, offset, h, us, UNIT_CONV_LENGTH);
+ writeArray(h_grp, fileName, xmfFile, "InternalEnergy", FLOAT, N, 1, parts,
+ N_total, mpi_rank, offset, u, us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N, 1, parts,
+ N_total, mpi_rank, offset, id, us, UNIT_CONV_NO_UNITS);
+ writeArray(h_grp, fileName, xmfFile, "Acceleration", FLOAT, N, 3, parts,
+ N_total, mpi_rank, offset, a_hydro, us, UNIT_CONV_ACCELERATION);
+ writeArray(h_grp, fileName, xmfFile, "Density", FLOAT, N, 1, parts, N_total,
+ mpi_rank, offset, rho, us, UNIT_CONV_DENSITY);
+}
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grpsph The HDF5 group in which to write
+ */
+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);
+
+ /* Viscosity and thermal conduction */
+ /* Nothing in this minimal model... */
+ writeAttribute_s(h_grpsph, "Thermal Conductivity Model", "No model");
+ 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, "Maximal Delta u change over dt",
+ const_max_u_change);
+}
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..b89cc249a42032e05bd35498c5c6435f3ebc4f4b
--- /dev/null
+++ b/src/hydro/Minimal/hydro_part.h
@@ -0,0 +1,108 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/**
+ * @brief Particle fields not needed during the SPH loops over neighbours.
+ *
+ * This structure contains the particle fields that are not used in the
+ * density or force loops. Quantities should be used in the kick, drift and
+ * potentially ghost tasks only.
+ */
+struct xpart {
+
+ double x_old[3]; /*!< Old position, at last tree rebuild. */
+
+ float v_full[3]; /*!< Velocity at the last full step. */
+
+ float u_full; /*!< Thermal energy at the last full step. */
+
+} __attribute__((aligned(xpart_align)));
+
+/**
+ * @brief Particle fields for the SPH particles
+ *
+ * The density and force substructures are used to contain variables only used
+ * within the density and force loops over neighbours. All more permanent
+ * variables should be declared in the main part of the part structure,
+ */
+struct part {
+
+ double x[3]; /*!< Particle position. */
+
+ float v[3]; /*!< Particle predicted velocity. */
+
+ float a_hydro[3]; /*!< Particle acceleration. */
+
+ float mass; /*!< Particle mass. */
+
+ float h; /*!< Particle smoothing length. */
+
+ float h_dt; /*!< Time derivative of smoothing length */
+
+ float t_begin; /*!< Time at the beginning of time-step. */
+
+ float t_end; /*!< Time at the end of time-step. */
+
+ float u; /*!< Particle internal energy. */
+
+ float u_dt; /*!< Time derivative of the internal energy. */
+
+ float rho; /*!< Particle density. */
+
+ float rho_dh; /*!< Derivative of density with respect to h */
+
+ /* Store density/force specific stuff. */
+ union {
+
+ /**
+ * @brief Structure for the variables only used in the density loop over
+ * neighbours.
+ *
+ * Quantities in this sub-structure should only be accessed in the density
+ * loop over neighbours and the ghost task.
+ */
+ struct {
+
+ float wcount; /*!< Neighbour number count. */
+
+ float wcount_dh; /*!< Derivative of the neighbour number with respect to
+ h. */
+ } density;
+
+ /**
+ * @brief Structure for the variables only used in the force loop over
+ * neighbours.
+ *
+ * Quantities in this sub-structure should only be accessed in the force
+ * loop over neighbours and the ghost and kick tasks.
+ */
+ struct {
+
+ float pressure; /*!< Particle pressure. */
+
+ float v_sig; /*!< Particle signal velocity */
+
+ } force;
+ };
+
+ unsigned long long id; /*!< Particle unique ID. */
+
+ struct gpart* gpart; /*!< Pointer to corresponding gravity part. */
+
+} __attribute__((aligned(part_align)));
diff --git a/src/hydro_io.h b/src/hydro_io.h
index b39661185668f71ae4adfe151de127a2a4766189..30d663f647c9b763e9b19177e9ba8ef374855768 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -22,10 +22,14 @@
#include "./const.h"
/* Import the right functions */
-#ifdef LEGACY_GADGET2_SPH
+#if defined(MINIMAL_SPH)
+#include "./hydro/Minimal/hydro_io.h"
+#elif defined(GADGET2_SPH)
#include "./hydro/Gadget2/hydro_io.h"
-#else
+#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_io.h"
+#else
+#error "Invalid choice of SPH variant"
#endif
-#endif
+#endif /* SWIFT_HYDRO_IO_H */
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 1391322af1e1520cabc007eee354b669bc4219c6..5ac0d7122b17c135dca386f3e9c74123d15676cd 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -459,7 +459,7 @@ void write_output_parallel(struct engine* e, struct UnitSystem* us,
int mpi_rank, int mpi_size, MPI_Comm comm,
MPI_Info info) {
- hid_t h_file = 0, h_grp = 0;
+ hid_t h_file = 0, h_grp = 0, h_grpsph = 0;
int N = e->s->nr_parts;
int periodic = e->s->periodic;
unsigned int numParticles[6] = {N, 0};
@@ -546,7 +546,10 @@ void write_output_parallel(struct engine* e, struct UnitSystem* us,
writeCodeDescription(h_file);
/* Print the SPH parameters */
- writeSPHflavour(h_file);
+ h_grpsph = H5Gcreate1(h_file, "/SPH", 0);
+ if (h_grpsph < 0) error("Error while creating SPH group");
+ writeSPHflavour(h_grpsph);
+ H5Gclose(h_grpsph);
/* Print the system of Units */
writeUnitSystem(h_file, us);
diff --git a/src/part.h b/src/part.h
index 3aaa28939c240bb24fb7053ad949f8b7ec999a88..168e80b68bb211d1d773f1f80f1fa9cc757edfcb 100644
--- a/src/part.h
+++ b/src/part.h
@@ -38,10 +38,14 @@
#define xpart_align 32
/* Import the right particle definition */
-#ifdef LEGACY_GADGET2_SPH
+#if defined(MINIMAL_SPH)
+#include "./hydro/Minimal/hydro_part.h"
+#elif defined(GADGET2_SPH)
#include "./hydro/Gadget2/hydro_part.h"
-#else
+#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_part.h"
+#else
+#error "Invalid choice of SPH variant"
#endif
#include "./gravity/Default/gravity_part.h"
diff --git a/src/runner.c b/src/runner.c
index bfdaba3c64d0a4589b6477f7c7fd98ebbb2797f2..157616ad29c954bb7ebfd72acfea96454327e4ee 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -34,6 +34,7 @@
#include "runner.h"
/* Local headers. */
+#include "approx_math.h"
#include "atomic.h"
#include "const.h"
#include "debug.h"
@@ -514,9 +515,6 @@ void runner_doinit(struct runner *r, struct cell *c, int timer) {
if (p->t_end <= t_end) {
- /* if(p->id == 1000) message("init 1000!"); */
- /* if(p->id == 515050) message("init 515050!"); */
-
/* Get ready for a density calculation */
hydro_init_part(p);
}
@@ -579,9 +577,6 @@ void runner_doghost(struct runner *r, struct cell *c) {
p = &parts[pid[i]];
xp = &xparts[pid[i]];
- /* if(p->id == 1000) message("ghost 1000"); */
- /* if(p->id == 515050) message("ghost 515050"); */
-
/* Is this part within the timestep? */
if (p->t_end <= t_end) {
@@ -619,14 +614,14 @@ void runner_doghost(struct runner *r, struct cell *c) {
}
/* We now have a particle whose smoothing length has converged */
- // if(p->id == 1000)
- // printParticle(parts, 1000, count);
- /* As of here, particle force variables will be set. Do _NOT_
- try to read any particle density variables! */
+ /* As of here, particle force variables will be set. */
/* Compute variables required for the force loop */
- hydro_prepare_force(p, xp);
+ hydro_prepare_force(p, xp, t_end);
+
+ /* The particle force values are now set. Do _NOT_
+ try to read any particle density variables! */
/* Prepare the particle for the force loop over neighbours */
hydro_reset_acceleration(p);
@@ -705,6 +700,7 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
struct part *restrict p, *restrict parts = c->parts;
struct xpart *restrict xp, *restrict xparts = c->xparts;
float dx_max = 0.f, h_max = 0.f;
+ float w;
TIMER_TIC
@@ -718,52 +714,36 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
p = &parts[k];
xp = &xparts[k];
+ /* Useful quantity */
+ const float h_inv = 1.0f / p->h;
+
/* Drift... */
p->x[0] += xp->v_full[0] * dt;
p->x[1] += xp->v_full[1] * dt;
p->x[2] += xp->v_full[2] * dt;
- /* Predict velocities */
- p->v[0] += p->a[0] * dt;
- p->v[1] += p->a[1] * dt;
- p->v[2] += p->a[2] * dt;
-
- /* /\* Predict smoothing length *\/ */
- /* w = p->force.h_dt * ih * dt; */
- /* if (fabsf(w) < 0.01f) /\* 1st order expansion of exp(w) *\/ */
- /* p->h *= */
- /* 1.0f + */
- /* w * (1.0f + w * (0.5f + w * (1.0f / 6.0f + 1.0f / 24.0f *
- * w))); */
- /* else */
- /* p->h *= expf(w); */
-
- // MATTHIEU
-
- /* /\* Predict density *\/ */
- /* w = -3.0f * p->force.h_dt * ih * dt; */
- /* if (fabsf(w) < 0.1f) */
- /* p->rho *= */
- /* 1.0f + */
- /* w * (1.0f + w * (0.5f + w * (1.0f / 6.0f + 1.0f / 24.0f *
- * w))); */
- /* else */
- /* p->rho *= expf(w); */
+ /* Predict velocities (for hydro terms) */
+ p->v[0] += p->a_hydro[0] * dt;
+ p->v[1] += p->a_hydro[1] * dt;
+ p->v[2] += p->a_hydro[2] * dt;
+
+ /* Predict smoothing length */
+ w = p->h_dt * h_inv * dt;
+ if (fabsf(w) < 0.2f)
+ p->h *= approx_expf(w); /* 4th order expansion of exp(w) */
+ else
+ p->h *= expf(w);
+
+ /* Predict density */
+ w = -3.0f * p->h_dt * h_inv * dt;
+ if (fabsf(w) < 0.2f)
+ p->rho *= approx_expf(w); /* 4th order expansion of exp(w) */
+ else
+ p->rho *= expf(w);
/* Predict the values of the extra fields */
hydro_predict_extra(p, xp, r->e->timeOld, r->e->time);
- /* if(p->id == 1000 || p->id == 515050 || p->id == 504849) */
- /* message("%lld: current_t=%f t0=%f t1=%f v=[%.3e %.3e %.3e]\n",
- */
- /* p->id, */
- /* r->e->time, */
- /* r->e->timeOld, */
- /* r->e->time, */
- /* p->v[0], */
- /* p->v[1], */
- /* p->v[2]); */
-
/* Compute motion since last cell construction */
const float dx =
sqrtf((p->x[0] - xp->x_old[0]) * (p->x[0] - xp->x_old[0]) +
@@ -827,9 +807,10 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
int updated = 0;
float t_end_min = FLT_MAX, t_end_max = 0.f;
- double ekin = 0.0, epot = 0.0;
- float mom[3] = {0.0f, 0.0f, 0.0f}, ang[3] = {0.0f, 0.0f, 0.0f};
- float m, x[3], v_full[3];
+ double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
+ float mom[3] = {0.0f, 0.0f, 0.0f};
+ float ang[3] = {0.0f, 0.0f, 0.0f};
+ float x[3], v_full[3];
struct part *restrict p, *restrict parts = c->parts;
struct xpart *restrict xp, *restrict xparts = c->xparts;
@@ -845,7 +826,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
p = &parts[k];
xp = &xparts[k];
- m = p->mass;
+ const float m = p->mass;
x[0] = p->x[0];
x[1] = p->x[1];
x[2] = p->x[2];
@@ -854,8 +835,13 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
if (is_fixdt || p->t_end <= t_current) {
/* First, finish the force loop */
+ p->h_dt *= p->h * 0.333333333f;
+
+ /* And do the same of the extra variable */
hydro_end_force(p);
+ /* Now we are ready to compute the next time-step size */
+
if (is_fixdt) {
/* Now we have a time step, proceed with the kick */
@@ -869,6 +855,13 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
new_dt = fminf(new_dt_hydro, new_dt_grav);
+ /* 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)
+ : FLT_MAX;
+
+ new_dt = fminf(new_dt, dt_h_change);
+
/* Recover the current timestep */
const float current_dt = p->t_end - p->t_begin;
@@ -898,28 +891,16 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
p->t_end = p->t_begin + new_dt;
/* Kick particles in momentum space */
- xp->v_full[0] += p->a[0] * dt;
- xp->v_full[1] += p->a[1] * dt;
- xp->v_full[2] += p->a[2] * dt;
-
- p->v[0] = xp->v_full[0] - half_dt * p->a[0];
- p->v[1] = xp->v_full[1] - half_dt * p->a[1];
- p->v[2] = xp->v_full[2] - half_dt * p->a[2];
-
- /* if(p->id == 1000 || p->id == 515050 || p->id == 504849) */
- /* message("%lld: current_t=%f t_beg=%f t_end=%f half_dt=%f v=[%.3e
- * %.3e %.3e]\n", */
- /* p->id, */
- /* t_current, */
- /* p->t_begin, */
- /* p->t_end, */
- /* half_dt, */
- /* p->v[0], */
- /* p->v[1], */
- /* p->v[2]); */
+ xp->v_full[0] += p->a_hydro[0] * dt;
+ xp->v_full[1] += p->a_hydro[1] * dt;
+ xp->v_full[2] += p->a_hydro[2] * dt;
+
+ p->v[0] = xp->v_full[0] - half_dt * p->a_hydro[0];
+ p->v[1] = xp->v_full[1] - half_dt * p->a_hydro[1];
+ p->v[2] = xp->v_full[2] - half_dt * p->a_hydro[2];
/* Extra kick work */
- hydro_kick_extra(p, dt);
+ hydro_kick_extra(p, xp, dt, half_dt);
}
/* Now collect quantities for statistics */
@@ -928,6 +909,9 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
v_full[1] = xp->v_full[1];
v_full[2] = xp->v_full[2];
+ /* Collect mass */
+ mass += m;
+
/* Collect momentum */
mom[0] += m * v_full[0];
mom[1] += m * v_full[1];
@@ -939,9 +923,10 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
ang[2] += m * (x[0] * v_full[1] - x[1] * v_full[0]);
/* Collect total energy. */
- ekin += 0.5 * m * (v_full[0] * v_full[0] + v_full[1] * v_full[1] +
- v_full[2] * v_full[2]);
- epot += m * xp->u_hdt;
+ e_kin += 0.5 * m * (v_full[0] * v_full[0] + v_full[1] * v_full[1] +
+ v_full[2] * v_full[2]);
+ e_pot += 0.f; /* No gravitational potential thus far */
+ e_int += hydro_get_internal_energy(p);
/* Minimal time for next end of time-step */
t_end_min = fminf(p->t_end, t_end_min);
@@ -960,11 +945,16 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
struct cell *cp = c->progeny[k];
+
+ /* Recurse */
runner_dokick(r, cp, 0);
+ /* And aggregate */
updated += cp->updated;
- ekin += cp->ekin;
- epot += cp->epot;
+ e_kin += cp->e_kin;
+ e_int += cp->e_int;
+ e_pot += cp->e_pot;
+ mass += cp->mass;
mom[0] += cp->mom[0];
mom[1] += cp->mom[1];
mom[2] += cp->mom[2];
@@ -978,8 +968,10 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* Store the values. */
c->updated = updated;
- c->ekin = ekin;
- c->epot = epot;
+ c->e_kin = e_kin;
+ c->e_int = e_int;
+ c->e_pot = e_pot;
+ c->mass = mass;
c->mom[0] = mom[0];
c->mom[1] = mom[1];
c->mom[2] = mom[2];
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index cddd249acd77f305a456d7b78023d4e4ca613cc7..1b987eb140569629fc4583199c58e3f33446cb78 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -1398,9 +1398,6 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
/* Get a pointer to the ith particle. */
pi = &parts[pid];
- /* if(pi->id == 1000) message("oO 1000"); */
- /* if(pi->id == 515050) message("oO 515050"); */
-
/* Get the particle position and radius. */
for (k = 0; k < 3; k++) pix[k] = pi->x[k];
hi = pi->h;
@@ -1624,9 +1621,6 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
/* Get a pointer to the ith particle. */
pi = &parts[pid];
- /* if(pi->id == 1000) message("oO 1000"); */
- /* if(pi->id == 515050) message("oO 515050"); */
-
/* Get the particle position and radius. */
for (k = 0; k < 3; k++) pix[k] = pi->x[k];
hi = pi->h;
diff --git a/src/serial_io.c b/src/serial_io.c
index ebb7eee68882766c5f550cda7feee1d9b68ab5cc..5402afba096bc62530ef04a5c717c2180d7a7c95 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -241,10 +241,10 @@ void prepareArray(hid_t grp, char* fileName, FILE* xmfFile, char* name,
*
* Calls #error() if an error occurs.
*/
-void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
- enum DATA_TYPE type, int N, int dim, long long N_total,
- int mpi_rank, long long offset, char* part_c,
- struct UnitSystem* us,
+void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
+ enum DATA_TYPE type, int N, int dim, long long N_total,
+ int mpi_rank, long long offset, char* part_c,
+ struct UnitSystem* us,
enum UnitConversionFactor convFactor) {
hid_t h_data = 0, h_err = 0, h_memspace = 0, h_filespace = 0;
@@ -259,11 +259,10 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
/* message("Writing '%s' array...", name); */
/* Prepare the arrays in the file */
- if(mpi_rank == 0)
- prepareArray(grp, fileName, xmfFile, name, type, N_total, dim,
- us, convFactor);
+ if (mpi_rank == 0)
+ prepareArray(grp, fileName, xmfFile, name, type, N_total, dim, us,
+ convFactor);
-
/* Allocate temporary buffer */
temp = malloc(N * dim * sizeOfType(type));
if (temp == NULL) error("Unable to allocate memory for temporary buffer");
@@ -359,9 +358,9 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
*/
#define writeArray(grp, fileName, xmfFile, name, type, N, dim, part, N_total, \
mpi_rank, offset, field, us, convFactor) \
- writeArrayBackEnd(grp, fileName, xmfFile, name, type, N, dim, N_total, \
- mpi_rank, offset, (char*)(&(part[0]).field), \
- us, convFactor)
+ writeArrayBackEnd(grp, fileName, xmfFile, name, type, N, dim, N_total, \
+ mpi_rank, offset, (char*)(&(part[0]).field), us, \
+ convFactor)
/* Import the right hydro definition */
#include "hydro_io.h"
@@ -512,7 +511,7 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, int* N,
*/
void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
int mpi_size, MPI_Comm comm, MPI_Info info) {
- hid_t h_file = 0, h_grp = 0;
+ hid_t h_file = 0, h_grp = 0, h_grpsph = 0;
int N = e->s->nr_parts;
int periodic = e->s->periodic;
int numParticles[6] = {N, 0};
@@ -602,7 +601,10 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
writeCodeDescription(h_file);
/* Print the SPH parameters */
- writeSPHflavour(h_file);
+ h_grpsph = H5Gcreate1(h_file, "/SPH", 0);
+ if (h_grpsph < 0) error("Error while creating SPH group");
+ writeSPHflavour(h_grpsph);
+ H5Gclose(h_grpsph);
/* Print the system of Units */
writeUnitSystem(h_file, us);
@@ -639,8 +641,8 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
error("Error while opening particle group on rank %d.\n", mpi_rank);
/* Write particle fields from the particle structure */
- hydro_write_particles(h_grp, fileName, xmfFile, N, N_total, mpi_rank, offset,
- parts, us);
+ hydro_write_particles(h_grp, fileName, xmfFile, N, N_total, mpi_rank,
+ offset, parts, us);
/* Close particle group */
H5Gclose(h_grp);
diff --git a/src/single_io.c b/src/single_io.c
index 1a74f690f1a0658f74a58c55d3ab0ffcf2717463..b0e350ffe2e2f7ff077c29dd0b581c9c3c555822 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -383,7 +383,7 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
*/
void write_output_single(struct engine* e, struct UnitSystem* us) {
- hid_t h_file = 0, h_grp = 0;
+ hid_t h_file = 0, h_grp = 0, h_grpsph = 0;
int N = e->s->nr_parts;
int periodic = e->s->periodic;
int numParticles[6] = {N, 0};
@@ -451,7 +451,10 @@ void write_output_single(struct engine* e, struct UnitSystem* us) {
writeCodeDescription(h_file);
/* Print the SPH parameters */
- writeSPHflavour(h_file);
+ h_grpsph = H5Gcreate1(h_file, "/SPH", 0);
+ if (h_grpsph < 0) error("Error while creating SPH group");
+ writeSPHflavour(h_grpsph);
+ H5Gclose(h_grpsph);
/* Print the system of Units */
writeUnitSystem(h_file, us);
diff --git a/src/tools.c b/src/tools.c
index 34cd4c436a128a3e177a7d4405b313130b143889..5feba7759f730faea1f38ceb9835f2076bc37a56 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -420,7 +420,7 @@ void engine_single_force(double *dim, long long int pid,
}
/* Dump the result. */
- message("part %lli (h=%e) has a=[%.3e,%.3e,%.3e]", p.id, p.h, p.a[0], p.a[1],
- p.a[2]);
+ message("part %lli (h=%e) has a=[%.3e,%.3e,%.3e]", p.id, p.h, p.a_hydro[0],
+ p.a_hydro[1], p.a_hydro[2]);
fflush(stdout);
}