diff --git a/examples/test_single.c b/examples/test_single.c
index f2aae8095d40d4ec5cb8e2617eb688b2039efa67..761b196df15051df38ebfc7e9c4fb67d64b6645c 100644
--- a/examples/test_single.c
+++ b/examples/test_single.c
@@ -78,10 +78,10 @@ int main ( int argc , char *argv[] ) {
p2.rho = 1.0f; p2.mass = 9.7059e-4; p2.h = 0.222871287 / 2;
p1.force.c = 0.0040824829f; p1.force.balsara = 0.0f;
p2.force.c = 58.8972740361f; p2.force.balsara = 0.0f;
- p1.u = 1.e-5 / ((const_gamma - 1.)*p1.rho);
- p2.u = 1.e-5 / ((const_gamma - 1.)*p2.rho) + 100.0f / ( 33 * p2.mass );
- p1.force.POrho2 = p1.u * ( const_gamma - 1.0f ) / p1.rho;
- p2.force.POrho2 = p2.u * ( const_gamma - 1.0f ) / p2.rho;
+ p1.u = 1.e-5 / ((const_hydro_gamma - 1.)*p1.rho);
+ p2.u = 1.e-5 / ((const_hydro_gamma - 1.)*p2.rho) + 100.0f / ( 33 * p2.mass );
+ p1.force.POrho2 = p1.u * ( const_hydro_gamma - 1.0f ) / p1.rho;
+ p2.force.POrho2 = p2.u * ( const_hydro_gamma - 1.0f ) / p2.rho;
/* Dump a header. */
printParticle_single( &p1 );
diff --git a/src/const.h b/src/const.h
index 09e0eb87a814f637cd4cc39f4f042e5615fc8ade..3caab0ba2b99b5396f36fd7ad0283a19deb51545 100644
--- a/src/const.h
+++ b/src/const.h
@@ -20,7 +20,7 @@
/* Hydrodynamical constants. */
-#define const_gamma (5.0f/3.0f)
+#define const_hydro_gamma (5.0f/3.0f)
#define const_viscosity_alpha 0.8f
/* Time integration constants. */
@@ -29,6 +29,6 @@
#define const_max_u_change 0.1f
/* Neighbour search constants. */
-#define const_eta_kernel 1.2348f /* Corresponds to 48 ngbs with the cubic spline kernel */
+#define const_eta_kernel 1.2349f /* Corresponds to 48 ngbs with the cubic spline kernel */
#define const_delta_nwneigh 1.f
#define CUBIC_SPLINE_KERNEL
diff --git a/src/engine.c b/src/engine.c
index 653255a3bca2ad3f4a1155232a732688832082fd..447c9c0709b9d9dedd6130bfb752a4e6894efd25 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -263,7 +263,7 @@ void engine_map_kick_first ( struct cell *c , void *data ) {
rho = p->rho *= 1.0f + w*( 1.0f + w*( -0.5f + w*(1.0f/6.0f - 1.0f/24.0*w ) ) );
else */
rho = p->rho *= expf( w );
- p->force.POrho2 = u * ( const_gamma - 1.0f ) / ( rho * rho );
+ p->force.POrho2 = u * ( const_hydro_gamma - 1.0f ) / ( rho * rho );
}
else {
p->density.wcount = 0.0f;
diff --git a/src/io.c b/src/io.c
index 4bde914a4a7d84b4271ded06e041ab584a4bb406..33a810850925fabb37500f8118322081cad52464 100644
--- a/src/io.c
+++ b/src/io.c
@@ -30,7 +30,7 @@
#include <string.h>
#include <stddef.h>
#include <hdf5.h>
-
+#include <math.h>
#include "lock.h"
#include "task.h"
@@ -38,13 +38,14 @@
#include "space.h"
#include "engine.h"
#include "error.h"
+#include "kernel.h"
/**
* @brief The different types of data used in the GADGET IC files.
*
* (This is admittedly a poor substitute to C++ templates...)
*/
-enum DATA_TYPE{INT, LONG, LONGLONG, UINT, ULONG, ULONGLONG, FLOAT, DOUBLE};
+enum DATA_TYPE{INT, LONG, LONGLONG, UINT, ULONG, ULONGLONG, FLOAT, DOUBLE, CHAR};
/**
* @brief The two sorts of data present in the GADGET IC files: compulsory to start a run or optional.
@@ -70,6 +71,7 @@ static hid_t hdf5Type(enum DATA_TYPE type)
case ULONGLONG: return H5T_NATIVE_ULLONG;
case FLOAT: return H5T_NATIVE_FLOAT;
case DOUBLE: return H5T_NATIVE_DOUBLE;
+ case CHAR: return H5T_C_S1;
default: error("Unknown type"); return 0;
}
}
@@ -89,6 +91,7 @@ static size_t sizeOfType(enum DATA_TYPE type)
case ULONGLONG: return sizeof(unsigned long long);
case FLOAT: return sizeof(float);
case DOUBLE: return sizeof(double);
+ case CHAR: return sizeof(char);
default: error("Unknown type"); return 0;
}
}
@@ -411,6 +414,120 @@ void writeAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data, int
H5Aclose(h_attr);
}
+/**
+ * @brief Write a string as an attribute to a given HDF5 group.
+ *
+ * @param grp The group in which to write.
+ * @param name The name of the attribute to write.
+ * @param str The string to write.
+ * @param length The length of the string
+ *
+ * Calls #error() if an error occurs.
+ */
+void writeStringAttribute(hid_t grp, char* name, char* str, int length)
+{
+ hid_t h_space=0, h_attr=0, h_err=0, h_type=0;
+
+ h_space = H5Screate(H5S_SCALAR);
+ if(h_space < 0)
+ {
+ char buf[100];
+ sprintf(buf, "Error while creating dataspace for attribute '%s'\n", name);
+ error(buf);
+ }
+
+ h_type = H5Tcopy(H5T_C_S1);
+ if(h_type < 0)
+ {
+ char buf[100];
+ sprintf(buf, "Error while copying datatype 'H5T_C_S1'\n");
+ error(buf);
+ }
+
+ h_err = H5Tset_size(h_type, length);
+ if(h_err < 0)
+ {
+ char buf[100];
+ sprintf(buf, "Error while resizing attribute tyep to '%i'\n", length);
+ error(buf);
+ }
+
+ h_attr = H5Acreate1(grp, name, h_type, h_space, H5P_DEFAULT);
+ if(h_attr < 0)
+ {
+ char buf[100];
+ sprintf(buf, "Error while creating attribute '%s'\n", name);
+ error(buf);
+ }
+
+ h_err = H5Awrite(h_attr, h_type, str );
+ if(h_err < 0)
+ {
+ char buf[100];
+ sprintf(buf, "Error while reading attribute '%s'\n", name);
+ error(buf);
+ }
+
+ H5Tclose(h_type);
+ H5Sclose(h_space);
+ H5Aclose(h_attr);
+}
+
+/**
+ * @brief Writes a double value as an attribute
+ * @param grp The groupm in which to write
+ * @param name The name of the attribute
+ * @param data The value to write
+ */
+void writeAttribute_d(hid_t grp, char* name, double data)
+{
+ writeAttribute(grp, name, DOUBLE, &data, 1);
+}
+
+/**
+ * @brief Writes a float value as an attribute
+ * @param grp The groupm in which to write
+ * @param name The name of the attribute
+ * @param data The value to write
+ */
+void writeAttribute_f(hid_t grp, char* name, float data)
+{
+ writeAttribute(grp, name, FLOAT, &data, 1);
+}
+
+/**
+ * @brief Writes an int value as an attribute
+ * @param grp The groupm in which to write
+ * @param name The name of the attribute
+ * @param data The value to write
+ */
+
+void writeAttribute_i(hid_t grp, char* name, int data)
+{
+ writeAttribute(grp, name, INT, &data, 1);
+}
+
+/**
+ * @brief Writes a long value as an attribute
+ * @param grp The groupm in which to write
+ * @param name The name of the attribute
+ * @param data The value to write
+ */
+void writeAttribute_l(hid_t grp, char* name, long data)
+{
+ writeAttribute(grp, name, LONG, &data, 1);
+}
+
+/**
+ * @brief Writes a string value as an attribute
+ * @param grp The groupm in which to write
+ * @param name The name of the attribute
+ * @param str The string to write
+ */
+void writeAttribute_s(hid_t grp, char* name, char* str)
+{
+ writeStringAttribute(grp, name, str, strlen(str));
+}
/**
@@ -542,7 +659,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, enu
void write_output (struct engine *e)
{
- 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};
@@ -588,8 +705,6 @@ void write_output (struct engine *e)
/* Close runtime parameters */
H5Gclose(h_grp);
-
- (void) sizeOfType(INT);
/* Open header to write simulation properties */
/* printf("write_output: Writing file header...\n"); */
@@ -597,7 +712,7 @@ void write_output (struct engine *e)
if(h_grp < 0)
error("Error while creating file header\n");
- /* Read the relevant information and print status */
+ /* Print the relevant information and print status */
writeAttribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 1);
writeAttribute(h_grp, "NumPart_ThisFile", UINT, numParticles, 6);
writeAttribute(h_grp, "Time", DOUBLE, &e->time, 1);
@@ -609,8 +724,25 @@ void write_output (struct engine *e)
writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, 6);
writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, numParticlesHighWord, 6);
writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
-
- /* Close header */
+
+ /* Print the SPH parameters */
+ h_grpsph = H5Gcreate1(h_file, "/Header/SPH", 0);
+ if(h_grpsph < 0)
+ error("Error while creating SPH header\n");
+
+ writeAttribute_f(h_grpsph, "Kernel eta", const_eta_kernel);
+ writeAttribute_f(h_grpsph, "Weighted N_ngb", kernel_nwneigh);
+ writeAttribute_f(h_grpsph, "Delta N_ngb", const_delta_nwneigh);
+ writeAttribute_f(h_grpsph, "Hydro gamma", const_hydro_gamma);
+ writeAttribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
+ 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);
+
+ /* Close headers */
+ H5Gclose(h_grpsph);
H5Gclose(h_grp);
/* Create SPH particles group */
diff --git a/src/kernel.h b/src/kernel.h
index 96912d2152c4e40e4747feeec7606da2bb7f2b26..35a204bca1e46bc346383442b28b299da8b8d6de 100644
--- a/src/kernel.h
+++ b/src/kernel.h
@@ -36,6 +36,7 @@
/* -------------------------------------------------------------------------------------------------------------------- */
/* Coefficients for the kernel. */
+#define kernel_name "Cubic spline"
#define kernel_degree 3
#define kernel_ivals 2
#define kernel_gamma 2.0f
@@ -127,6 +128,7 @@ __attribute__ ((always_inline)) INLINE static void kernel_eval ( float x , float
/* -------------------------------------------------------------------------------------------------------------------- */
/* Coefficients for the kernel. */
+#define kernel_name "Quartic spline"
#define kernel_degree 4
#define kernel_ivals 3
#define kernel_gamma 2.5f
@@ -220,6 +222,7 @@ __attribute__ ((always_inline)) INLINE static void kernel_eval ( float x , float
/* -------------------------------------------------------------------------------------------------------------------- */
/* Coefficients for the kernel. */
+#define kernel_name "Quintic spline"
#define kernel_degree 5
#define kernel_ivals 3
#define kernel_gamma 3.f
diff --git a/src/runner.c b/src/runner.c
index d8c032b3917a30f5a528bd9c959f6cf72b84c175..785bb72e2248105f9811ba2928d983a287048085 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -426,13 +426,13 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
/* Compute this particle's sound speed. */
u = p->u;
- p->force.c = fc = sqrtf( const_gamma * ( const_gamma - 1.0f ) * u );
+ p->force.c = fc = sqrtf( const_hydro_gamma * ( const_hydro_gamma - 1.0f ) * u );
/* Compute the P/Omega/rho2. */
- p->force.POrho2 = u * ( const_gamma - 1.0f ) / ( rho + h * rho_dh / 3.0f );
+ p->force.POrho2 = u * ( const_hydro_gamma - 1.0f ) / ( rho + h * rho_dh / 3.0f );
- /* Balsara switch */
- p->force.balsara = normDiv_v / ( normDiv_v + normCurl_v + 0.0001f * fc * ih );
+ /* Balsara switch */
+ p->force.balsara = normDiv_v / ( normDiv_v + normCurl_v + 0.0001f * fc * ih );
/* Reset the acceleration. */
for ( k = 0 ; k < 3 ; k++ )