From 79e36b809fa186c02f8d56a7195c9d7f2f46b5f7 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Fri, 30 Nov 2012 16:13:13 +0000
Subject: [PATCH] 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
---
examples/Makefile.am | 2 +
examples/test.c | 254 ++++++++++++++++++++++++++----------------
src/Makefile.am | 9 +-
src/gadgetsmp.h | 2 +-
src/ic.c | 260 ++++++++++++++++++++++++++++++++++++++++++-
src/ic.h | 7 +-
6 files changed, 435 insertions(+), 99 deletions(-)
diff --git a/examples/Makefile.am b/examples/Makefile.am
index ce5007d377..361de521c0 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -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
diff --git a/examples/test.c b/examples/test.c
index 0498a9d024..5bcec2a09b 100644
--- a/examples/test.c
+++ b/examples/test.c
@@ -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 );
diff --git a/src/Makefile.am b/src/Makefile.am
index addcde017e..33907eb666 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -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
diff --git a/src/gadgetsmp.h b/src/gadgetsmp.h
index 3193cc8200..67a01137a8 100644
--- a/src/gadgetsmp.h
+++ b/src/gadgetsmp.h
@@ -32,4 +32,4 @@
#include "runner_iact.h"
#include "engine.h"
#include "runner.h"
-
+#include "ic.h"
diff --git a/src/ic.c b/src/ic.c
index 315395fa26..11025277e2 100644
--- a/src/ic.c
+++ b/src/ic.c
@@ -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 */
+
diff --git a/src/ic.h b/src/ic.h
index c604b9a66b..10216bfb8c 100644
--- a/src/ic.h
+++ b/src/ic.h
@@ -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);
+
+
--
GitLab