Commit 79e36b80 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

The ICs for particles and the box can now be read from a (GADGET-)HDF5 file....

The ICs for particles and the box can now be read from a (GADGET-)HDF5 file. The file name is passed to the comand line via the option -f and all the (now) unnecessary comand line options have been removed. All the values that can be read from the files are read from the file. The source file ic.c can be easily changed to read additional arrays from the HDF5 inputs. The Makefiles have been updated to reflect the changes and the need to link against the (serial) HDF5 libreary.


Former-commit-id: 5b8896489301cbad20fe566c72e624b52a8f593f
parent 6b7df171
......@@ -21,6 +21,8 @@ AUTOMAKE_OPTIONS=gnu
# Add the source directory and debug to CFLAGS
AM_CFLAGS = -g -Wall -Werror -I../src $(OPENMP_CFLAGS) -DCPU_TPS=2.67e9
AM_LDFLAGS = $(HDF5_LDFLAGS)
# Set-up the library
bin_PROGRAMS = test
......
......@@ -628,108 +628,176 @@ int main ( int argc , char *argv[] ) {
int nr_threads = 1, nr_queues = -1, runs = 1;
int data[2];
double dim[3] = { 1.0 , 1.0 , 1.0 }, shift[3] = { 0.0 , 0.0 , 0.0 };
double r_min = 0.01, r_max = 0.1, h_max = -1.0 , scaling = 1.0, rho = 0.0;
double /*r_min = 0.01, r_max = 0.1,*/ h_max = -1.0 , scaling = 1.0, rho = 0.0;
struct part *parts = NULL, *p;
struct space s;
struct engine e;
char ICfileName[200];
ticks tic;
/* Init the space. */
bzero( &s , sizeof(struct space) );
/* Parse the options. */
while ( ( c = getopt( argc , argv , "a:b:p:d:N:c:h:v:m:s:t:q:r:i:m:z:" ) ) != -1 )
switch ( c ) {
case 'N':
if ( sscanf( optarg , "%d" , &N ) != 1 )
error( "Error parsing number of particles." );
if ( posix_memalign( (void *)&parts , 32 , N * sizeof(struct part) ) != 0 )
error( "Call to posix_memalign failed." );
for ( k = 0 ; k < N ; k++ ) {
parts[k].x[0] = ((double)rand()) / RAND_MAX * dim[0];
parts[k].x[1] = ((double)rand()) / RAND_MAX * dim[1];
parts[k].x[2] = ((double)rand()) / RAND_MAX * dim[2];
parts[k].id = k;
parts[k].h = r_min + ((r_max - r_min)*rand())/RAND_MAX;
parts[k].mass = 1.0f;
parts[k].u = 1.0f;
}
printf( "main: allocated memory for %i parts.\n" , N ); fflush(stdout);
break;
case 'a':
if ( sscanf( optarg , "%lf" , &scaling ) != 1 )
error( "Error parsing cutoff scaling." );
printf( "main: scaling cutoff by %.3f.\n" , scaling ); fflush(stdout);
for ( k = 0 ; k < N ; k++ )
parts[k].h *= scaling;
break;
case 'b':
if ( sscanf( optarg , "%lf %lf %lf" , &dim[0] , &dim[1] , &dim[2] ) != 3 )
error( "Error parsing box dimensions." );
break;
case 'c':
printf( "main: reading parts from %s...\n" , optarg ); fflush(stdout);
if ( parts == NULL && posix_memalign( (void *)&parts , 16 , N * sizeof(struct part) ) != 0 )
error( "Call to calloc failed." );
read_coords( optarg , parts , N );
break;
case 'd':
printf( "main: reading dt from %s...\n" , optarg ); fflush(stdout);
read_dt( optarg , parts , N );
break;
case 'h':
printf( "main: reading cutoffs from %s...\n" , optarg ); fflush(stdout);
read_cutoffs( optarg , parts , N );
break;
case 'i':
printf( "main: reading ids from %s...\n" , optarg ); fflush(stdout);
read_id( optarg , parts , N );
break;
case 'm':
if ( sscanf( optarg , "%lf" , &h_max ) != 1 )
error( "Error parsing h_max." );
printf( "main: maximum h set to %e.\n" , h_max );
break;
case 'p':
if ( sscanf( optarg , "%d" , &periodic ) != 1 )
error( "Error parsing periodicity." );
printf( "main: periodicity switched %s.\n" , periodic ? "on" : "off" );
break;
case 'q':
if ( sscanf( optarg , "%d" , &nr_queues ) != 1 )
error( "Error parsing number of queues." );
break;
case 'r':
if ( sscanf( optarg , "%d" , &runs ) != 1 )
error( "Error parsing number of runs." );
break;
case 's':
if ( sscanf( optarg , "%lf %lf %lf" , &shift[0] , &shift[1] , &shift[2] ) != 3 )
error( "Error parsing shift." );
for ( k = 0 ; k < N ; k++ ) {
parts[k].x[0] += shift[0];
parts[k].x[1] += shift[1];
parts[k].x[2] += shift[2];
}
printf( "main: shifted parts by [ %.3f %.3f %.3f ].\n" , shift[0] , shift[1] , shift[2] );
break;
case 't':
if ( sscanf( optarg , "%d" , &nr_threads ) != 1 )
error( "Error parsing number of threads." );
omp_set_num_threads( nr_threads );
break;
case 'z':
if ( sscanf( optarg , "%d" , &space_splitsize ) != 1 )
error( "Error parsing split size." );
printf( "main: split size set to %i.\n" , space_splitsize );
break;
case '?':
error( "Unknown option." );
break;
}
/* Parse the options */
while ( ( c = getopt( argc , argv , "f:a:m:s:t:q:r:z:" ) ) != -1 )
switch( c )
{
case 'f':
if( !strcpy(ICfileName, optarg))
error("Error parsing IC file name.");
printf("main: IC to be read from file '%s'\n", ICfileName);
break;
case 'a':
if ( sscanf( optarg , "%lf" , &scaling ) != 1 )
error( "Error parsing cutoff scaling." );
printf( "main: scaling cutoff by %.3f.\n" , scaling ); fflush(stdout);
break;
case 'm':
if ( sscanf( optarg , "%lf" , &h_max ) != 1 )
error( "Error parsing h_max." );
printf( "main: maximum h set to %e.\n" , h_max ); fflush(stdout);
break;
case 'q':
if ( sscanf( optarg , "%d" , &nr_queues ) != 1 )
error( "Error parsing number of queues." );
break;
case 'r':
if ( sscanf( optarg , "%d" , &runs ) != 1 )
error( "Error parsing number of runs." );
break;
case 's':
if ( sscanf( optarg , "%lf %lf %lf" , &shift[0] , &shift[1] , &shift[2] ) != 3 )
error( "Error parsing shift." );
printf( "main: will shift parts by [ %.3f %.3f %.3f ].\n" , shift[0] , shift[1] , shift[2] );
break;
case 't':
if ( sscanf( optarg , "%d" , &nr_threads ) != 1 )
error( "Error parsing number of threads." );
omp_set_num_threads( nr_threads );
break;
case 'z':
if ( sscanf( optarg , "%d" , &space_splitsize ) != 1 )
error( "Error parsing split size." );
printf( "main: split size set to %i.\n" , space_splitsize );
break;
case '?':
error( "Unknown option." );
break;
}
/* while ( ( c = getopt( argc , argv , "a:b:p:d:N:c:h:v:m:s:t:q:r:i:m:z:" ) ) != -1 ) */
/* switch ( c ) { */
/* case 'N': */
/* if ( sscanf( optarg , "%d" , &N ) != 1 ) */
/* error( "Error parsing number of particles." ); */
/* if ( posix_memalign( (void *)&parts , 32 , N * sizeof(struct part) ) != 0 ) */
/* error( "Call to posix_memalign failed." ); */
/* for ( k = 0 ; k < N ; k++ ) { */
/* parts[k].x[0] = ((double)rand()) / RAND_MAX * dim[0]; */
/* parts[k].x[1] = ((double)rand()) / RAND_MAX * dim[1]; */
/* parts[k].x[2] = ((double)rand()) / RAND_MAX * dim[2]; */
/* parts[k].id = k; */
/* parts[k].h = r_min + ((r_max - r_min)*rand())/RAND_MAX; */
/* parts[k].mass = 1.0f; */
/* parts[k].u = 1.0f; */
/* } */
/* printf( "main: allocated memory for %i parts.\n" , N ); fflush(stdout); */
/* break; */
/* case 'a': */
/* if ( sscanf( optarg , "%lf" , &scaling ) != 1 ) */
/* error( "Error parsing cutoff scaling." ); */
/* printf( "main: scaling cutoff by %.3f.\n" , scaling ); fflush(stdout); */
/* for ( k = 0 ; k < N ; k++ ) */
/* parts[k].h *= scaling; */
/* break; */
/* case 'b': */
/* if ( sscanf( optarg , "%lf %lf %lf" , &dim[0] , &dim[1] , &dim[2] ) != 3 ) */
/* error( "Error parsing box dimensions." ); */
/* break; */
/* case 'c': */
/* printf( "main: reading parts from %s...\n" , optarg ); fflush(stdout); */
/* if ( parts == NULL && posix_memalign( (void *)&parts , 16 , N * sizeof(struct part) ) != 0 ) */
/* error( "Call to calloc failed." ); */
/* read_coords( optarg , parts , N ); */
/* break; */
/* case 'd': */
/* printf( "main: reading dt from %s...\n" , optarg ); fflush(stdout); */
/* read_dt( optarg , parts , N ); */
/* break; */
/* case 'h': */
/* printf( "main: reading cutoffs from %s...\n" , optarg ); fflush(stdout); */
/* read_cutoffs( optarg , parts , N ); */
/* break; */
/* case 'i': */
/* printf( "main: reading ids from %s...\n" , optarg ); fflush(stdout); */
/* read_id( optarg , parts , N ); */
/* break; */
/* case 'm': */
/* if ( sscanf( optarg , "%lf" , &h_max ) != 1 ) */
/* error( "Error parsing h_max." ); */
/* printf( "main: maximum h set to %e.\n" , h_max ); */
/* break; */
/* case 'p': */
/* if ( sscanf( optarg , "%d" , &periodic ) != 1 ) */
/* error( "Error parsing periodicity." ); */
/* printf( "main: periodicity switched %s.\n" , periodic ? "on" : "off" ); */
/* break; */
/* case 'q': */
/* if ( sscanf( optarg , "%d" , &nr_queues ) != 1 ) */
/* error( "Error parsing number of queues." ); */
/* break; */
/* case 'r': */
/* if ( sscanf( optarg , "%d" , &runs ) != 1 ) */
/* error( "Error parsing number of runs." ); */
/* break; */
/* case 's': */
/* if ( sscanf( optarg , "%lf %lf %lf" , &shift[0] , &shift[1] , &shift[2] ) != 3 ) */
/* error( "Error parsing shift." ); */
/* for ( k = 0 ; k < N ; k++ ) { */
/* parts[k].x[0] += shift[0]; */
/* parts[k].x[1] += shift[1]; */
/* parts[k].x[2] += shift[2]; */
/* } */
/* printf( "main: shifted parts by [ %.3f %.3f %.3f ].\n" , shift[0] , shift[1] , shift[2] ); */
/* break; */
/* case 't': */
/* if ( sscanf( optarg , "%d" , &nr_threads ) != 1 ) */
/* error( "Error parsing number of threads." ); */
/* omp_set_num_threads( nr_threads ); */
/* break; */
/* case 'z': */
/* if ( sscanf( optarg , "%d" , &space_splitsize ) != 1 ) */
/* error( "Error parsing split size." ); */
/* printf( "main: split size set to %i.\n" , space_splitsize ); */
/* break; */
/* case '?': */
/* error( "Unknown option." ); */
/* break; */
/* } */
/* How large are the parts? */
printf( "main: sizeof(struct part) is %li bytes.\n" , (long int)sizeof( struct part ) );
printf( "main: sizeof(struct part) is %li bytes.\n" , (long int)sizeof( struct part ));
/* Read particles and space information from (GADGET) IC */
tic = getticks();
read_ic(ICfileName, dim, &parts, &N, &periodic);
printf( "main: reading particle properties took %.3f ms.\n" , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);
/* Apply h scaling */
if(scaling != 1.0)
for ( k = 0 ; k < N ; k++ )
parts[k].h *= scaling;
/* Apply shift */
if(shift[0] !=0 || shift[1] !=0 || shift[2] !=0 )
for ( k = 0 ; k < N ; k++ ) {
parts[k].x[0] += shift[0];
parts[k].x[1] += shift[1];
parts[k].x[2] += shift[2];
}
/* Dump the kernel to make sure its ok. */
// kernel_dump( 100 );
......
......@@ -22,7 +22,7 @@ AUTOMAKE_OPTIONS=gnu
AM_CFLAGS = -g -O3 -Wall -Werror -ffast-math -fstrict-aliasing $(SIMD_FLAGS) $(CFLAGS) $(OPENMP_CFLAGS) -DTIMER -DCOUNTER -DCPU_TPS=2.67e9
# Assign a "safe" version number
AM_LDFLAGS = $(LAPACK_LIBS) $(BLAS_LIBS) -version-info 0:0:0
AM_LDFLAGS = $(LAPACK_LIBS) $(BLAS_LIBS) $(HDF5_LDFLAGS) -version-info 0:0:0
# Build the libgadgetsmp library
lib_LTLIBRARIES = libgadgetsmp.la
......@@ -32,3 +32,10 @@ libgadgetsmp_la_SOURCES = space.c runner.c queue.c task.c cell.c engine.c ic.c
include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
engine.h gadgetsmp.h ic.h
# Set-up the library
#bin_PROGRAMS = test
# Sources for test
#test_SOURCES = main.c
#test_CFLAGS = -DCOUNTER -DTIMER $(AM_CFLAGS)
#test_LDADD = ../src/.libs/libgadgetsmp.a
......@@ -32,4 +32,4 @@
#include "runner_iact.h"
#include "engine.h"
#include "runner.h"
#include "ic.h"
......@@ -22,31 +22,287 @@
#ifdef HAVE_HDF5
/* Some standard headers. */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stddef.h>
#include <hdf5.h>
#include "task.h"
#include "lock.h"
#include "part.h"
#include "space.h"
/* Error macro. */
#define error(s) { printf( "%s:%s():%i: %s\n" , __FILE__ , __FUNCTION__ , __LINE__ , s ); abort(); }
/**
* @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};
/**
* @brief Converts a C data type to the HDF5 equivalent.
*
* This function is a trivial wrapper around the HDF5 types but allows
* to change the exact storage types matching the code types in a transparent way.
*/
hid_t hdf5Type(enum DATA_TYPE type)
{
switch(type)
{
case INT: return H5T_NATIVE_INT;
case UINT: return H5T_NATIVE_UINT;
case LONG: return H5T_NATIVE_LONG;
case ULONG: return H5T_NATIVE_ULONG;
case LONGLONG: return H5T_NATIVE_LLONG;
case ULONGLONG: return H5T_NATIVE_ULLONG;
case FLOAT: return H5T_NATIVE_FLOAT;
case DOUBLE: return H5T_NATIVE_DOUBLE;
default: error("Unknown type");
}
}
/**
* @brief Returns the memory size of the data type
*/
size_t sizeOfType(enum DATA_TYPE type)
{
switch(type)
{
case INT: return sizeof(int);
case UINT: return sizeof(unsigned int);
case LONG: return sizeof(long);
case ULONG: return sizeof(unsigned long);
case LONGLONG: return sizeof(long long);
case ULONGLONG: return sizeof(unsigned long long);
case FLOAT: return sizeof(float);
case DOUBLE: return sizeof(double);
default: error("Unknown type");
}
}
/**
* @brief Reads an attribute from a given HDF5 group.
*
* @param grp The group from which to read.
* @param name The name of the attribute to read.
* @param type The #DATA_TYPE of the attribute.
* @param data (output) The attribute read from the HDF5 group.
*
* Calls #error() if an error occurs.
*/
void readAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data)
{
hid_t h_attr=0, h_err=0;
h_attr = H5Aopen(grp, name, H5P_DEFAULT);
if(h_attr < 0)
{
char buf[100];
sprintf(buf, "Error while opening attribute '%s'\n", name);
error(buf);
}
h_err = H5Aread(h_attr, hdf5Type(type), data);
if(h_err < 0)
{
char buf[100];
sprintf(buf, "Error while reading attribute '%s'\n", name);
error(buf);
}
H5Aclose(h_attr);
}
/**
* @brief Reads a data array from a given HDF5 group.
*
* @param grp The group from which to read.
* @param name The name of the array to read.
* @param type The #DATA_TYPE of the attribute.
* @param N The number of particles.
* @param dim The dimension of the data (1 for scalar, 3 for vector)
* @param part_c A (char*) pointer on the first occurence of the field of interest in the parts array
*
* @todo A better version using HDF5 hyperslabs to read the file directly into the part array
* will be written once the strucutres have been stabilized.
*
* Calls #error() if an error occurs.
*/
void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N, int dim, char* part_c)
{
hid_t h_data=0, h_err=0, h_type=0;
void* temp;
int i=0;
size_t typeSize = sizeOfType(type);
size_t copySize = typeSize * dim;
size_t partSize = sizeof(struct part);
char* temp_c = 0;
printf("readArray: Reading '%s' array...\n", name);
/* Open data space */
h_data = H5Dopen(grp, name);
if(h_data < 0)
{
char buf[100];
sprintf(buf, "Error while opening data space '%s'\n", name);
error(buf);
}
/* Check data type */
h_type = H5Dget_type(h_data);
if(h_type < 0)
error("Unable to retrieve data type from the file");
if(!H5Tequal(h_type, hdf5Type(type)))
error("Non-matching types between the code and the file");
/* Allocate temporary buffer */
temp = malloc(N * dim * sizeOfType(type));
if(temp == NULL)
error("Unable to allocate memory for temporary buffer");
/* Read HDF5 dataspace in temporary buffer */
/* Dirty version that happens to work for vectors but should be improved */
/* Using HDF5 dataspaces would be better */
h_err = H5Dread(h_data, hdf5Type(type), H5S_ALL, H5S_ALL, H5P_DEFAULT, temp);
if(h_err < 0)
{
char buf[100];
sprintf(buf, "Error while reading data array '%s'\n", name);
error(buf);
}
/* Copy temporary buffer to particle data */
temp_c = temp;
for(i=0; i<N; ++i)
memcpy(part_c+i*partSize, &temp_c[i*copySize], copySize);
/* Free and close everything */
free(temp);
H5Tclose(h_type);
H5Dclose(h_data);
}
/**
* @brief A helper macro to call the readArrayBackEnd function more easily.
*
* @param grp The group from which to read.
* @param name The name of the array to read.
* @param type The #DATA_TYPE of the attribute.
* @param N The number of particles.
* @param dim The dimension of the data (1 for scalar, 3 for vector)
* @param part The array of particles to fill
* @param field The name of the field (C code name as defined in part.h) to fill
*
*/
#define readArray(grp, name, type, N, dim, part, field) readArrayBackEnd(grp, name, type, N, dim, (char*)(&(part[0]).field))
/**
* @brief Reads an HDF5 initial condition file
* @brief Reads an HDF5 initial condition file (GADGET-3 type)
*
* @param fileName The file to read.
* @param dim (output) The dimension of the volume read from the file.
* @param parts (output) The array of #part read from the file.
* @param N (output) The number of particles read from the file.
* @param periodic (output) 1 if the volume is periodic, 0 if not.
*
* Opens the HDF5 file fileName and reads the particles contained
* in the parts array. N is the returned number of particles found
* in the file.
*
* @warning Can not read snapshot distributed over more than 1 file !!!
* @todo Read snaphsots distributed in more than one file.
*
* Calls #error() if an error occurs.
*
*/
void read_ic ( char* fileName, struct part *parts, int* N)
void read_ic ( char* fileName, double dim[3], struct part **parts, int* N, int* periodic)
{
hid_t h_file=0, h_grp=0;
double boxSize=0.; /* GADGET has only cubic boxes (in cosmological mode) */
int numParticles[6]={0}; /* GADGET has 6 particle types. We only keep the type 0*/
/* Open file */
printf("read_ic: Opening file '%s' as IC.\n", fileName);
h_file = H5Fopen(fileName, H5F_ACC_RDONLY, H5P_DEFAULT);
if(h_file < 0)
{
char buf[200];
sprintf(buf, "Error while opening file '%s'", fileName);
error(buf);
}
/* Open header to read simulation properties */
printf("read_ic: Reading runtime parameters...\n");
h_grp = H5Gopen(h_file, "/RuntimePars");
if(h_grp < 0)
error("Error while opening runtime parameters\n");
/* Read the relevant information */
readAttribute(h_grp, "PeriodicBoundariesOn", INT, periodic);
/* Close runtime parameters */
H5Gclose(h_grp);
/* Open header to read simulation properties */
printf("read_ic: Reading file header...\n");
h_grp = H5Gopen(h_file, "/Header");
if(h_grp < 0)
error("Error while opening file header\n");
/* Read the relevant information and print status */
readAttribute(h_grp, "BoxSize", DOUBLE, &boxSize);
readAttribute(h_grp, "NumPart_Total", UINT, numParticles);
*N = numParticles[0];
dim[0] = dim[1] = dim[2] = boxSize;
printf("read_ic: Found %d particles in a %speriodic box of size [%f %f %f]\n",
*N, (periodic ? "": "non-"), dim[0], dim[1], dim[2]);
/* Close header */
H5Gclose(h_grp);
/* Allocate memory to store particles */
if(posix_memalign( (void*)parts , 32 , *N * sizeof(struct part)) != 0)
error("Error while allocating memory for particles");
printf("read_ic: Allocated %8.2f MB for particles.\n", *N * sizeof(struct part) / (1024.*1024.));
/* Open SPH particles group */
printf("read_ic: Reading particle arrays...\n");
h_grp = H5Gopen(h_file, "/PartType0");
if(h_grp < 0)
error( "Error while opening particle group.\n");
/* Read arrays */
readArray(h_grp, "Coordinates", DOUBLE, *N, 3, *parts, x);
readArray(h_grp, "Velocity", FLOAT, *N, 3, *parts, v);
readArray(h_grp, "Mass", FLOAT, *N, 1, *parts, mass);
readArray(h_grp, "SmoothingLength", FLOAT, *N, 1, *parts, h);
readArray(h_grp, "InternalEnergy", FLOAT, *N, 1, *parts, u);
readArray(h_grp, "ParticleIDs", ULONG, *N, 1, *parts, id);
/* Close particle group */
H5Gclose(h_grp);
printf("read_ic: Done Reading particles...\n");
/* Close file */
H5Fclose(h_file);
}
#endif /* HAVE_HDF5 */
......@@ -18,5 +18,8 @@
******************************************************************************/
void read_ic ( char* fileName, struct space *s, int* N);
void space_init ( struct space *s , double dim[3] , struct part *parts , int N , int periodic , double h_max );
void read_ic ( char* fileName, double dim[3], struct part **parts, int* N, int* periodic);
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