Commit 9907bd7f authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

The different constants used in SPH are now written to the HDF5 output for...

The different constants used in SPH are now written to the HDF5 output for better tracability. The const_gamma constant has been renamed to const_hydro_gamma.


Former-commit-id: 2f117589f48ad150e47d62733c4d3635caf74295
parent 628978c1
......@@ -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 );
......
......@@ -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
......@@ -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;
......
......@@ -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 */
......
......@@ -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
......
......@@ -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++ )
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment