io.c 28.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
/*******************************************************************************
 * This file is part of SWIFT.
 * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
 *                    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/>.
 * 
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

#ifdef HAVE_HDF5


/* Some standard headers. */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stddef.h>
#include <hdf5.h>
33
#include <math.h>
34
35
36
37
38

#include "lock.h"
#include "task.h"
#include "part.h"
#include "space.h"
39
#include "engine.h"
40
#include "error.h"
41
#include "kernel.h"
42
43
44
45
46
47

/**
 * @brief The different types of data used in the GADGET IC files.
 *
 * (This is admittedly a poor substitute to C++ templates...)
 */
48
enum DATA_TYPE{INT, LONG, LONGLONG, UINT, ULONG, ULONGLONG, FLOAT, DOUBLE, CHAR};
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

/**
 * @brief The two sorts of data present in the GADGET IC files: compulsory to start a run or optional.
 *
 */
enum DATA_IMPORTANCE{COMPULSORY=1, OPTIONAL=0};

/**
 * @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.
 */
static 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;
74
    case CHAR: return H5T_C_S1;
75
    default: error("Unknown type"); return 0;
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
    }
}

/**
 * @brief Returns the memory size of the data type
 */
static 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);
94
    case CHAR: return sizeof(char);
95
    default: error("Unknown type"); return 0;
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
    }
}

/*-----------------------------------------------------------------------------
 * Routines reading an IC file
 *-----------------------------------------------------------------------------*/


/**
 * @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
 * @param importance If COMPULSORY, the data must be present in the IC file. If OPTIONAL, the array will be zeroed when the data is not present.
 *
 * @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, enum DATA_IMPORTANCE importance)
{
  hid_t h_data=0, h_err=0, h_type=0;
  htri_t exist=0;
  void* temp;
  int i=0;
  const size_t typeSize = sizeOfType(type);
  const size_t copySize = typeSize * dim;
  const size_t partSize = sizeof(struct part);
  char* temp_c = 0;

  /* Check whether the dataspace exists or not */
  exist = H5Lexists(grp, name, 0);
  if(exist < 0)
    {
      char buf[100];
      sprintf(buf, "Error while checking the existence of data set '%s'\n", name);
      error(buf);
    }
  else if(exist == 0)
    {
      if(importance == COMPULSORY)
	{
	  char buf[100];
	  sprintf(buf, "Compulsory data set '%s' not present in the file.\n", name);
	  error(buf);
	}
      else
	{
	  /* printf("readArray: Optional data set '%s' not present. Zeroing this particle field...\n", name);	   */
	  
	  for(i=0; i<N; ++i)
	    memset(part_c+i*partSize, 0, copySize);
	  
	  return;
	}
   }

  /* printf("readArray: Reading %s '%s' array...\n", importance == COMPULSORY ? "compulsory": "optional  ", name); */

  /* Open data space */
  h_data = H5Dopen1(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
 * @param importance Is the data compulsory or not
 *
 */
#define readArray(grp, name, type, N, dim, part, field, importance) readArrayBackEnd(grp, name, type, N, dim, (char*)(&(part[0]).field), importance)

/**
 * @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, double dim[3], struct part **parts,  int* N, int* periodic)
{
  hid_t h_file=0, h_grp=0;
274
  double boxSize[3]={0.0,-1.0,-1.0};         /* GADGET has only cubic boxes (in cosmological mode) */
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
  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 = H5Gopen1(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 = H5Gopen1(h_file, "/Header");
  if(h_grp < 0)
    error("Error while opening file header\n");
    
  /* Read the relevant information and print status */
Pedro Gonnet's avatar
Pedro Gonnet committed
306
  readAttribute(h_grp, "BoxSize", DOUBLE, boxSize);
307
308
309
  readAttribute(h_grp, "NumPart_Total", UINT, numParticles);

  *N = numParticles[0];
Pedro Gonnet's avatar
Pedro Gonnet committed
310
  dim[0] = boxSize[0];
311
312
  dim[1] = ( boxSize[1] < 0 ) ? boxSize[0] : boxSize[1];
  dim[2] = ( boxSize[2] < 0 ) ? boxSize[0] : boxSize[2];
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334

  /* 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");
  bzero( *parts , *N * sizeof(struct part) );

  /* 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 = H5Gopen1(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, COMPULSORY);
335
336
  readArray(h_grp, "Velocities", FLOAT, *N, 3, *parts, v, COMPULSORY);
  readArray(h_grp, "Masses", FLOAT, *N, 1, *parts, mass, COMPULSORY);
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
  readArray(h_grp, "SmoothingLength", FLOAT, *N, 1, *parts, h, COMPULSORY);
  readArray(h_grp, "InternalEnergy", FLOAT, *N, 1, *parts, u, COMPULSORY);
  readArray(h_grp, "ParticleIDs", ULONGLONG, *N, 1, *parts, id, COMPULSORY);
  readArray(h_grp, "TimeStep", FLOAT, *N, 1, *parts, dt, OPTIONAL);
  readArray(h_grp, "Acceleration", FLOAT, *N, 3, *parts, a, OPTIONAL);
  readArray(h_grp, "Density", FLOAT, *N, 1, *parts, rho, OPTIONAL );

  /* Close particle group */
  H5Gclose(h_grp);

  /* printf("read_ic: Done Reading particles...\n"); */

  /* Close file */
  H5Fclose(h_file);
}


/*-----------------------------------------------------------------------------
 * Routines writing an output file
 *-----------------------------------------------------------------------------*/
357
358


359
/* File descriptor output routine */
360
361
void createXMFfile();
FILE* prepareXMFfile();
362
void writeXMFfooter(FILE* xmfFile);
363
void writeXMFheader(FILE* xmfFile, int Nparts, char* hdfFileName,  float time);
364
365
366
void writeXMFline(FILE* xmfFile, char* fileName, char* name, int N, int dim, enum DATA_TYPE type);


367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
/**
 * @brief Write 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 type The #DATA_TYPE of the attribute.
 * @param data The attribute to write.
 * @param num The number of elements to write
 *
 * Calls #error() if an error occurs.
 */
void writeAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data, int num)
{
  hid_t h_space=0, h_attr=0, h_err=0;
  hsize_t dim[1]={num};

  h_space = H5Screate(H5S_SIMPLE);
  if(h_space < 0)
    {
      char buf[100];
      sprintf(buf, "Error while creating dataspace for attribute '%s'\n", name);
      error(buf);
    }

  h_err = H5Sset_extent_simple(h_space, 1, dim, NULL);
  if(h_err < 0)
    {
      char buf[100];
      sprintf(buf, "Error while changing dataspace shape for attribute '%s'\n", name);
      error(buf);
    }

399
  h_attr = H5Acreate1(grp, name, hdf5Type(type), h_space, H5P_DEFAULT);
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
  if(h_attr < 0)
    {
      char buf[100];
      sprintf(buf, "Error while creating attribute '%s'\n", name);
      error(buf);
    }

  h_err = H5Awrite(h_attr, hdf5Type(type), data);
  if(h_err < 0)
    {
      char buf[100];
      sprintf(buf, "Error while reading attribute '%s'\n", name);
      error(buf);
    }

  H5Sclose(h_space);
  H5Aclose(h_attr);
}

419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
/**
 * @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));
}
533
534
535
536
537
538


/**
 * @brief Writes a data array in given HDF5 group.
 *
 * @param grp The group in which to write.
539
540
 * @param fileName The name of the file in which the data is written
 * @param xmfFile The FILE used to write the XMF description
541
542
543
544
545
546
547
548
549
550
551
 * @param name The name of the array to write.
 * @param type The #DATA_TYPE of the array.
 * @param N The number of particles to write.
 * @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 write the file directly from the part array
 * will be written once the strucutres have been stabilized.
 *
 * Calls #error() if an error occurs.
 */
552
void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, enum DATA_TYPE type, int N, int dim, char* part_c)
553
554
555
556
557
558
559
560
561
562
{
  hid_t h_data=0, h_err=0, h_space=0;
  void* temp = 0;
  int i=0, rank=0;
  const size_t typeSize = sizeOfType(type);
  const size_t copySize = typeSize * dim;
  const size_t partSize = sizeof(struct part);
  char* temp_c = 0;
  hsize_t shape[2];

563
  /* printf("writeArray: Writing '%s' array...\n", name); */
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604

  /* Allocate temporary buffer */
  temp = malloc(N * dim * sizeOfType(type));
  if(temp == NULL)
    error("Unable to allocate memory for temporary buffer");

  /* Copy temporary buffer to particle data */
  temp_c = temp;
  for(i=0; i<N; ++i)
    memcpy(&temp_c[i*copySize], part_c+i*partSize, copySize);

  /* Create data space */
  h_space = H5Screate(H5S_SIMPLE);
  if(h_space < 0)
    {
      char buf[100];
      sprintf(buf, "Error while creating data space for field '%s'\n", name);
      error(buf);      
    }
  
  if(dim > 1)
    {
      rank = 2;
      shape[0] = N; shape[1] = dim;
    }
  else
    {
      rank = 1;
      shape[0] = N; shape[1] = 0;
    }
  
  /* Change shape of data space */
  h_err = H5Sset_extent_simple(h_space, rank, shape, NULL);
  if(h_err < 0)
    {
      char buf[100];
      sprintf(buf, "Error while changing data space shape for field '%s'\n", name);
      error(buf);      
    }
  
  /* Create dataset */
605
  h_data = H5Dcreate1(grp, name, hdf5Type(type), h_space, H5P_DEFAULT);
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
  if(h_data < 0)
    {
      char buf[100];
      sprintf(buf, "Error while creating dataspace '%s'\n", name);
      error(buf);
    }
  
  /* Write temporary buffer to HDF5 dataspace */
  h_err = H5Dwrite(h_data, hdf5Type(type), h_space, H5S_ALL, H5P_DEFAULT, temp);
  if(h_err < 0)
    {
      char buf[100];
      sprintf(buf, "Error while reading data array '%s'\n", name);
      error(buf);
    }
621
622
623
624

  /* Write XMF description for this data set */
  writeXMFline(xmfFile, fileName, name, N, dim, type);

625
626
627
628
629
630
631
632
633
634
635
  
  /* Free and close everything */
  free(temp);
  H5Dclose(h_data);
  H5Sclose(h_space);
}

/**
 * @brief A helper macro to call the readArrayBackEnd function more easily.
 *
 * @param grp The group in which to write.
636
637
 * @param fileName The name of the file in which the data is written
 * @param xmfFile The FILE used to write the XMF description
638
639
640
641
 * @param name The name of the array to write.
 * @param type The #DATA_TYPE of the array.
 * @param N The number of particles to write.
 * @param dim The dimension of the data (1 for scalar, 3 for vector)
642
643
 * @param part A (char*) pointer on the first occurence of the field of interest in the parts array
 * @param field The name (code name) of the field to read from.
644
645
 *
 */
646
#define writeArray(grp, fileName, xmfFile, name, type, N, dim, part, field) writeArrayBackEnd(grp, fileName, xmfFile, name, type, N, dim, (char*)(&(part[0]).field))
647
648

/**
649
 * @brief Writes an HDF5 output file (GADGET-3 type) with its XMF descriptor
650
 *
651
 * @param e The engine containing all the system.
652
 *
653
 * Creates an HDF5 output file and writes the particles contained
654
 * in the engine. If such a file already exists, it is erased and replaced
655
656
 * by the new one. 
 * The companion XMF file is also updated accordingly.
657
658
659
660
 *
 * Calls #error() if an error occurs.
 *
 */
661
void write_output (struct engine *e)
662
{
663
  
664
  hid_t h_file=0, h_grp=0, h_grpsph=0;
665
666
  int N = e->s->nr_parts;
  int periodic = e->s->periodic;
667
  int numParticles[6]={N,0};
668
669
670
  int numParticlesHighWord[6]={0};
  int numFiles = 1;
  struct part* parts = e->s->parts;
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
  FILE* xmfFile = 0;
  static int outputCount = 0;
  
  /* File name */
  char fileName[200];
  sprintf(fileName, "output_%03i.hdf5", outputCount);

  /* First time, we need to create the XMF file */
  if(outputCount == 0)
    createXMFfile();
  
  /* Prepare the XMF file for the new entry */
  xmfFile = prepareXMFfile();

  /* Write the part corresponding to this specific output */
  writeXMFheader(xmfFile, N, fileName, e->time);
687

688
689

  /* Open file */
690
  /* printf("write_output: Opening file '%s'.\n", fileName); */
691
692
693
694
695
696
697
698
699
  h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT,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 */
700
  /* printf("write_output: Writing runtime parameters...\n"); */
701
  h_grp = H5Gcreate1(h_file, "/RuntimePars", 0);
702
703
704
705
706
707
708
709
710
711
  if(h_grp < 0)
    error("Error while creating runtime parameters group\n");

  /* Write the relevant information */
  writeAttribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);

  /* Close runtime parameters */
  H5Gclose(h_grp);
  
  /* Open header to write simulation properties */
712
  /* printf("write_output: Writing file header...\n"); */
713
  h_grp = H5Gcreate1(h_file, "/Header", 0);
714
715
716
  if(h_grp < 0)
    error("Error while creating file header\n");
    
717
  /* Print the relevant information and print status */
718
719
720
  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);
721

722
723
724
  /* GADGET-2 legacy values */
  writeAttribute(h_grp, "NumPart_Total", UINT, numParticles, 6);
  writeAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticlesHighWord, 6);
725
726
  double MassTable[6] = {0., 0., 0., 0., 0., 0.};
  writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, 6);
727
728
  writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, numParticlesHighWord, 6);
  writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747

  /* 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);
748
749
750
  H5Gclose(h_grp);
		  
  /* Create SPH particles group */
751
  /* printf("write_output: Writing particle arrays...\n"); */
752
  h_grp = H5Gcreate1(h_file, "/PartType0", 0);
753
754
755
756
  if(h_grp < 0)
    error( "Error while creating particle group.\n");

  /* Write arrays */
757
  writeArray(h_grp, fileName, xmfFile, "Coordinates", DOUBLE, N, 3, parts, x);
758
759
  writeArray(h_grp, fileName, xmfFile, "Velocities", FLOAT, N, 3, parts, v);
  writeArray(h_grp, fileName, xmfFile, "Masses", FLOAT, N, 1, parts, mass);
760
761
762
763
764
765
  writeArray(h_grp, fileName, xmfFile, "SmoothingLength", FLOAT, N, 1, parts, h);
  writeArray(h_grp, fileName, xmfFile, "InternalEnergy", FLOAT, N, 1, parts, u);
  writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N, 1, parts, id);
  writeArray(h_grp, fileName, xmfFile, "TimeStep", FLOAT, N, 1, parts, dt);
  writeArray(h_grp, fileName, xmfFile, "Acceleration", FLOAT, N, 3, parts, a);
  writeArray(h_grp, fileName, xmfFile, "Density", FLOAT, N, 1, parts, rho);
766
767
768
769

  /* Close particle group */
  H5Gclose(h_grp);

770
771
772
  /* Write LXMF file descriptor */
  writeXMFfooter(xmfFile);

773
  /* printf("write_output: Done writing particles...\n"); */
774
775
776

  /* Close file */
  H5Fclose(h_file);
777
778

  ++outputCount;
779
780
}

781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808


/* ------------------------------------------------------------------------------------------------ 
 * This part writes the XMF file descriptor enabling a visualisation through ParaView
 * ------------------------------------------------------------------------------------------------ */

/**
 * @brief Writes the lines corresponding to an array of the HDF5 output
 *
 * @param xmfFile The file in which to write
 * @param fileName The name of the HDF5 file associated to this XMF descriptor.
 * @param name The name of the array in the HDF5 file.
 * @param N The number of particles.
 * @param dim The dimension of the quantity (1 for scalars, 3 for vectors).
 * @param type The type of the data to write.
 *
 * @todo Treat the types in a better way.
 */
void writeXMFline(FILE* xmfFile, char* fileName, char* name, int N, int dim, enum DATA_TYPE type )
{
  fprintf(xmfFile, "<Attribute Name=\"%s\" AttributeType=\"%s\" Center=\"Node\">\n", name, dim == 1 ? "Scalar": "Vector");
  if(dim == 1)
    fprintf(xmfFile, "<DataItem Dimensions=\"%d\" NumberType=\"Double\" Precision=\"%d\" Format=\"HDF\">%s:/PartType0/%s</DataItem>\n", N, type==FLOAT ? 4:8, fileName, name);
  else
    fprintf(xmfFile, "<DataItem Dimensions=\"%d %d\" NumberType=\"Double\" Precision=\"%d\" Format=\"HDF\">%s:/PartType0/%s</DataItem>\n", N, dim, type==FLOAT ? 4:8, fileName, name);
  fprintf(xmfFile, "</Attribute>\n");
}

809

810
/**
811
 * @brief Prepares the XMF file for the new entry
812
 *
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
 * Creates a temporary file on the disk in order to copy the right lines.
 *
 * @todo Use a proper XML library to avoid stupid copies.
 */
FILE* prepareXMFfile()
{
  char buffer[1024];

  FILE* xmfFile = fopen("output.xmf", "r");
  FILE* tempFile = fopen("output_temp.xmf", "w");

  if(xmfFile == NULL)
    error("Unable to open current XMF file.");

  if(tempFile == NULL)
    error("Unable to open temporary file.");


  /* First we make a temporary copy of the XMF file and count the lines */
  int counter = 0;
  while (fgets(buffer, 1024, xmfFile) != NULL)
    {
      counter++;
      fprintf(tempFile, "%s", buffer);
    }
  fclose(tempFile);
  fclose(xmfFile);
  
  /* We then copy the XMF file back up to the closing lines */
  xmfFile = fopen("output.xmf", "w");
  tempFile = fopen("output_temp.xmf", "r");

  if(xmfFile == NULL)
    error("Unable to open current XMF file.");

  if(tempFile == NULL)
    error("Unable to open temporary file.");

  int i = 0;
  while (fgets(buffer, 1024, tempFile) != NULL && i < counter - 3)
    {
      i++;
      fprintf(xmfFile, "%s", buffer);
    }
  fprintf(xmfFile, "\n");
  fclose(tempFile);
  remove("output_temp.xmf");
 
  return xmfFile;
}

/**
 * @brief Writes the begin of the XMF file
866
867
868
 *
 * @todo Exploit the XML nature of the XMF format to write a proper XML writer and simplify all the XMF-related stuff.
 */
869
void createXMFfile()
870
{
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
  FILE* xmfFile = fopen("output.xmf", "w");

  fprintf(xmfFile, "<?xml version=\"1.0\" ?> \n");
  fprintf(xmfFile, "<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []> \n");
  fprintf(xmfFile, "<Xdmf xmlns:xi=\"http://www.w3.org/2003/XInclude\" Version=\"2.1\">\n");
  fprintf(xmfFile, "<Domain>\n");
  fprintf(xmfFile, "<Grid Name=\"TimeSeries\" GridType=\"Collection\" CollectionType=\"Temporal\">\n\n");

  fprintf(xmfFile, "</Grid>\n");
  fprintf(xmfFile, "</Domain>\n");
  fprintf(xmfFile, "</Xdmf>\n");

  fclose(xmfFile);
}

886

887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
/**
 * @brief Writes the part of the XMF entry presenting the geometry of the snapshot
 *
 * @param xmfFile The file to write in.
 * @param Nparts The number of particles.
 * @param hdfFileName The name of the HDF5 file corresponding to this output.
 * @param time The current simulation time.
 */
void writeXMFheader(FILE* xmfFile, int Nparts, char* hdfFileName, float time)
{
  /* Write end of file */
  
  fprintf(xmfFile, "<Grid GridType=\"Collection\" CollectionType=\"Spatial\">\n");
  fprintf(xmfFile, "<Time Type=\"Single\" Value=\"%f\"/>\n", time);
  fprintf(xmfFile, "<Grid Name=\"Gas\" GridType=\"Uniform\">\n");
  fprintf(xmfFile, "<Topology TopologyType=\"Polyvertex\" Dimensions=\"%d\"/>\n", Nparts);
  fprintf(xmfFile, "<Geometry GeometryType=\"XYZ\">\n");
  fprintf(xmfFile, "<DataItem Dimensions=\"%d 3\" NumberType=\"Double\" Precision=\"8\" Format=\"HDF\">%s:/PartType0/Coordinates</DataItem>\n", Nparts, hdfFileName);
  fprintf(xmfFile, "</Geometry>");
906
907
}

908

909
910
911
912
913
914
915
/**
 * @brief Writes the end of the XMF file (closes all open markups)
 *
 * @param xmfFile The file to write in.
 */
void writeXMFfooter(FILE* xmfFile)
{
916
  /* Write end of the section of this time step */
917
  
918
  fprintf(xmfFile, "\n</Grid>\n");
919
  fprintf(xmfFile, "</Grid>\n");
920
  fprintf(xmfFile, "\n</Grid>\n");
921
922
923
924
925
926
  fprintf(xmfFile, "</Domain>\n");
  fprintf(xmfFile, "</Xdmf>\n");
  
  fclose(xmfFile);
}

927
928
929
#endif  /* HAVE_HDF5 */