read_write_scalar_field.c 12.2 KB
Newer Older
gsell's avatar
gsell committed
1
/*
gsell's avatar
gsell committed
2
  Copyright (c) 2006-2017, The Regents of the University of California,
gsell's avatar
gsell committed
3 4 5 6 7 8
  through Lawrence Berkeley National Laboratory (subject to receipt of any
  required approvals from the U.S. Dept. of Energy) and the Paul Scherrer
  Institut (Switzerland).  All rights reserved.

  License: see file COPYING in top level of source distribution.
*/
9 10

#include "H5hut.h"
gsell's avatar
gsell committed
11
#include "examples.h"
12

13 14
#include <stdlib.h>

15
#define NPROCS  8
gsell's avatar
gsell committed
16
#define DEFAULT_VERBOSITY       H5_VERBOSE_DEFAULT
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 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 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 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

struct H5BlockPartition {
	h5_int64_t	i_start;
	h5_int64_t	i_end;
	h5_int64_t	j_start;
	h5_int64_t	j_end;
	h5_int64_t	k_start;
	h5_int64_t	k_end;
};

struct H5BlockPartition Layout1[1] = {
	{ 0, 63, 0, 63, 0, 511 }
};

struct H5BlockPartition Layout8[8] = {
		{  0,63,  0,63,   0, 63},
		{  0,63,  0,63,  64,127},
		{  0,63,  0,63, 128,191},
		{  0,63,  0,63, 192,255},
		{  0,63,  0,63, 256,319},
		{  0,63,  0,63, 320,383},
		{  0,63,  0,63, 384,447},
		{  0,63,  0,63, 448,511}
};

struct H5BlockPartition Layout8G[8] = {
		{  0,63,  0,63,   0, 64},
		{  0,63,  0,63,  63,128},
		{  0,63,  0,63, 127,192},
		{  0,63,  0,63, 191,256},
		{  0,63,  0,63, 255,320},
		{  0,63,  0,63, 319,384},
		{  0,63,  0,63, 383,448},
		{  0,63,  0,63, 447,511}
};

struct H5BlockPartition Layout16[16] = {
		{  0,63,  0,31,   0, 63},
		{  0,63, 32,63,   0, 63},
		{  0,63,  0,31,  64,127},
		{  0,63, 32,63,  64,127},
		{  0,63,  0,31, 128,191},
		{  0,63, 32,63, 128,191},
		{  0,63,  0,31, 192,255},
		{  0,63, 32,63, 192,255},
		{  0,63,  0,31, 256,319},
		{  0,63, 32,63, 256,319},
		{  0,63,  0,31, 320,383},
		{  0,63, 32,63, 320,383},
		{  0,63,  0,31, 384,447},
		{  0,63, 32,63, 384,447},
		{  0,63,  0,31, 448,511},
		{  0,63, 32,63, 448,511}
};

struct H5BlockPartition Layout16G[16] = {
		{  0,63,  0,32,   0, 64},
		{  0,63, 31,63,   0, 64},
		{  0,63,  0,32,  63,128},
		{  0,63, 31,63,  63,128},
		{  0,63,  0,32, 127,192},
		{  0,63, 31,63, 127,192},
		{  0,63,  0,32, 191,256},
		{  0,63, 31,63, 191,256},
		{  0,63,  0,32, 255,320},
		{  0,63, 31,63, 255,320},
		{  0,63,  0,32, 319,384},
		{  0,63, 31,63, 319,384},
		{  0,63,  0,32, 383,448},
		{  0,63, 31,63, 383,448},
		{  0,63,  0,32, 447,511},
		{  0,63, 31,63, 447,511}
};


struct H5BlockPartition Layout32[32] = {
		{  0,31,  0,31,   0, 63},
		{  0,31, 32,63,   0, 63},
		{ 32,63,  0,31,   0, 63},
		{ 32,63, 32,63,   0, 63},
		{  0,31,  0,31,  64,127},
		{  0,31, 32,63,  64,127},
		{ 32,63,  0,31,  64,127},
		{ 32,63, 32,63,  64,127},
		{  0,31,  0,31, 128,191},
		{  0,31, 32,63, 128,191},
		{ 32,63,  0,31, 128,191},
		{ 32,63, 32,63, 128,191},
		{  0,31,  0,31, 192,255},
		{  0,31, 32,63, 192,255},
		{ 32,63,  0,31, 192,255},
		{ 32,63, 32,63, 192,255},
		{  0,31,  0,31, 256,319},
		{  0,31, 32,63, 256,319},
		{ 32,63,  0,31, 256,319},
		{ 32,63, 32,63, 256,319},
		{  0,31,  0,31, 320,383},
		{  0,31, 32,63, 320,383},
		{ 32,63,  0,31, 320,383},
		{ 32,63, 32,63, 320,383},
		{  0,31,  0,31, 384,447},
		{  0,31, 32,63, 384,447},
		{ 32,63,  0,31, 384,447},
		{ 32,63, 32,63, 384,447},
		{  0,31,  0,31, 448,511},
		{  0,31, 32,63, 448,511},
		{ 32,63,  0,31, 448,511},
		{ 32,63, 32,63, 448,511}
};

struct H5BlockPartition Layout32G[32] = {
		{  0,32,  0,32,   0, 64},
		{  0,32, 31,63,   0, 64},
		{ 31,63,  0,32,   0, 64},
		{ 31,63, 31,63,   0, 64},
		{  0,32,  0,32,  63,128},
		{  0,32, 31,63,  63,128},
		{ 31,63,  0,32,  63,128},
		{ 31,63, 31,63,  63,128},
		{  0,32,  0,32, 127,192},
		{  0,32, 31,63, 127,192},
		{ 31,63,  0,32, 127,192},
		{ 31,63, 31,63, 127,192},
		{  0,32,  0,32, 191,256},
		{  0,32, 31,63, 191,256},
		{ 31,63,  0,32, 191,256},
		{ 31,63, 31,63, 191,256},
		{  0,32,  0,32, 255,320},
		{  0,32, 31,63, 255,320},
		{ 31,63,  0,32, 255,320},
		{ 31,63, 31,63, 255,320},
		{  0,32,  0,32, 319,384},
		{  0,32, 31,63, 319,384},
		{ 31,63,  0,32, 319,384},
		{ 31,63, 31,63, 319,384},
		{  0,31,  0,31, 383,448},
		{  0,31, 31,63, 383,448},
		{ 31,63,  0,31, 383,448},
		{ 31,63, 31,63, 383,448},
		{  0,32,  0,32, 447,511},
		{  0,32, 31,63, 447,511},
		{ 31,63,  0,32, 447,511},
		{ 31,63, 31,63, 447,511}
};

#define _calc_index( i, i_dims, j, j_dims, k, k_dims ) \
		(i + j*i_dims + k*i_dims*j_dims)

static h5_int64_t
_write_data (
	h5_file_t f,
	int myproc,
        struct H5BlockPartition* view
	) {

	h5_int64_t i, j, k, idx;
	h5_int64_t herr;
	h5_float64_t *data;
	h5_int64_t i_dims = view->i_end - view->i_start + 1;
	h5_int64_t j_dims = view->j_end - view->j_start + 1;
	h5_int64_t k_dims = view->k_end - view->k_start + 1;

	printf ( "Writing scalar field data to step #%lld\n", (long long)H5GetStep (f));

181
	data = (h5_float64_t*)malloc ( i_dims * j_dims * k_dims * sizeof (h5_float64_t) );
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
	for ( i = 0; i < i_dims; i++ ) {
		for ( j = 0; j < j_dims; j++ ) {
			for ( k = 0; k < k_dims; k++ ) {
				idx = _calc_index (
					i, i_dims,
					j, j_dims,
					k, k_dims );
				*(data + idx) = k
					+ 1000*j
					+ 100000*i
					+ 10000000*myproc;
			}
		}
	}

	herr = H5Block3dSetView (
		f,
		view->i_start, view->i_end,
		view->j_start, view->j_end,
		view->k_start, view->k_end );
	if ( herr < 0 ) return herr;

	herr = H5Block3dWriteScalarFieldFloat64 ( f, "TestField", data );
	if ( herr < 0 ) return herr;

	free ( data );

	return 1;
}

static h5_int64_t
_write_attributes (
	h5_file_t f,
	const int myproc
	) {
217 218
        printf ("Writing attributes to field '%s' in step #%lld\n",
		"TestField", (long long)H5GetStep (f));
219 220 221 222 223 224 225
	h5_int64_t herr = H5BlockWriteFieldAttribString (
		f,
		"TestField",
		"TestString",
		"42" );
	if ( herr < 0 ) return -1;

226 227 228 229 230 231 232 233 234
	h5_int32_t i4_val[1] = { 42 };
	herr = H5BlockWriteFieldAttribInt32 (
		f,
		"TestField",
		"TestInt32",
		i4_val, 1 );
	if ( herr < 0 ) return -1;

	h5_int64_t i8_val[1] = { 42 };
235 236 237 238
	herr = H5BlockWriteFieldAttribInt64 (
		f,
		"TestField",
		"TestInt64",
239 240 241 242 243 244 245 246 247
		i8_val, 1 );
	if ( herr < 0 ) return -1;

	h5_float32_t r4_val[1] = { 42.0 };
	herr = H5BlockWriteFieldAttribFloat32 (
		f,
		"TestField",
		"TestFloat32",
		r4_val, 1 );
248 249
	if ( herr < 0 ) return -1;

250
	h5_float64_t r8_val[1] = { 42.0 };
251 252 253 254
	herr = H5BlockWriteFieldAttribFloat64 (
		f,
		"TestField",
		"TestFloat64",
255
		r8_val, 1 );
256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273
	if ( herr < 0 ) return -1;

	herr = H5Block3dSetFieldOrigin ( f, "TestField", 1.0, 2.0, 3.0 );
	if ( herr < 0 ) return -1;

	herr = H5Block3dSetFieldSpacing ( f, "TestField", 2.0, 3.0, 4.0 );
	if ( herr < 0 ) return -1;

	return H5_SUCCESS;
}

static h5_int64_t
_write_file (
	const char *fname,
	const int myproc,
	MPI_Comm comm,
	struct H5BlockPartition *layout
	) {
gsell's avatar
gsell committed
274
		printf ("PROC[%d]: Open file \"%s\" for writing ...\n",
275
		myproc, fname );
gsell's avatar
gsell committed
276 277
        h5_file_t file = H5OpenFile (fname, H5_O_WRONLY, H5_PROP_DEFAULT);
        H5SetStep (file, 0);
278
	
gsell's avatar
gsell committed
279 280
	_write_data (file, myproc, layout);
	_write_attributes (file, myproc);
281
	
gsell's avatar
gsell committed
282
	H5CloseFile (file);
283
	
gsell's avatar
gsell committed
284
	return H5_SUCCESS;
285 286 287 288 289 290 291 292 293 294 295 296 297 298 299
}

static h5_int64_t
_read_data (
	h5_file_t f,
	int myproc,
	struct H5BlockPartition *layout
	) {

	h5_int64_t i, j, k, idx;
	h5_float64_t *data;
	h5_int64_t i_dims = layout->i_end - layout->i_start + 1;
	h5_int64_t j_dims = layout->j_end - layout->j_start + 1;
	h5_int64_t k_dims = layout->k_end - layout->k_start + 1;

gsell's avatar
gsell committed
300
	printf ("Reading Step #%lld\n", (long long)H5GetStep (f));
301

302
	data = (h5_float64_t*)malloc (i_dims * j_dims * k_dims * sizeof (*data));
303

gsell's avatar
gsell committed
304
	H5Block3dSetView (
305 306 307 308
		f,
		layout->i_start, layout->i_end,
		layout->j_start, layout->j_end,
		layout->k_start, layout->k_end );
gsell's avatar
gsell committed
309
	H5Block3dReadScalarFieldFloat64 ( f, "TestField", data );
310

gsell's avatar
gsell committed
311 312 313
	for (i = 0; i < i_dims; i++) {
		for (j = 0; j < j_dims; j++) {
			for (k = 0; k < k_dims; k++) {
314 315 316 317 318 319 320 321 322
				idx = _calc_index (
					i, i_dims,
					j, j_dims,
					k, k_dims );

				h5_float64_t value = k
					+ 1000 * j
					+ 100000 * i
					+ 10000000 * myproc;
gsell's avatar
gsell committed
323
				if (*(data + idx) != value) {
324 325 326 327 328 329 330 331 332 333 334 335 336
					printf (
						"PROC[%d]: "
						"value missmatch for (%lld,%lld,%lld); is: %f;"
						" should be: %f\n",
						myproc,
						(long long)i, (long long)j, (long long)k,
						*( data + idx ), value );
					return -1;
				}
			}
		}
	}

gsell's avatar
gsell committed
337
	free (data);
338

gsell's avatar
gsell committed
339
	return H5_SUCCESS;
340 341 342 343 344 345 346 347 348 349
}

static h5_int64_t
_read_attributes (
	h5_file_t f,
	const int myproc,
	MPI_Comm comm
	) {
	h5_int64_t timestep = 0;
	
gsell's avatar
gsell committed
350
	H5SetStep (f, timestep);
351 352

	char sval[16];
gsell's avatar
gsell committed
353
	H5BlockReadFieldAttribString (
354 355 356 357
		f,
		"TestField",
		"TestString",
		sval );
gsell's avatar
gsell committed
358 359 360
	if (strcmp (sval, "42") != 0) {
		printf ("Error reading string attribute: "
			"Value is \"%s\" and should be \"42\"\n", sval);
361 362
	}

363 364 365 366 367 368 369 370 371 372 373 374 375
	h5_int32_t i4_val[1];
	H5BlockReadFieldAttribInt32 (
		f,
		"TestField",
		"TestInt32",
		i4_val );
	if (i4_val[0] != 42) {
		printf ("Error reading int32 attribute: "
			"Value is %lld and should be 42\n",
			(long long) i4_val[0]);
	}

	h5_int64_t i8_val[1];
gsell's avatar
gsell committed
376
	H5BlockReadFieldAttribInt64 (
377 378 379
		f,
		"TestField",
		"TestInt64",
380 381
		i8_val );
	if (i8_val[0] != 42) {
gsell's avatar
gsell committed
382 383
		printf ("Error reading int64 attribute: "
			"Value is %lld and should be 42\n",
384 385 386 387 388 389 390 391 392 393 394 395 396
			(long long) i8_val[0]);
	}

	h5_float32_t r4_val[1];
	H5BlockReadFieldAttribFloat32 (
		f,
		"TestField",
		"TestFloat32",
		r4_val );
	if (r4_val[0] != 42.0) {
		printf ("Error reading float64 attribute: "
			"Value is %f and should be 42.0\n",
			r4_val[0]);
397 398
	}

399
	h5_float64_t r8_val[1];
gsell's avatar
gsell committed
400
	H5BlockReadFieldAttribFloat64 (
401 402 403
		f,
		"TestField",
		"TestFloat64",
404 405
		r8_val );
	if (r8_val[0] != 42.0) {
gsell's avatar
gsell committed
406 407
		printf ("Error reading float64 attribute: "
			"Value is %f and should be 42.0\n",
408
			r8_val[0]);
409 410 411 412 413 414 415 416 417
	}

	h5_float64_t x_origin;
	h5_float64_t y_origin;
	h5_float64_t z_origin;
	h5_float64_t x_spacing;
	h5_float64_t y_spacing;
	h5_float64_t z_spacing;

gsell's avatar
gsell committed
418
	H5Block3dGetFieldOrigin (
419 420 421 422 423
		f, "TestField",
		&x_origin,
		&y_origin,
		&z_origin );

gsell's avatar
gsell committed
424
	if (x_origin != 1.0 || y_origin != 2.0 || z_origin != 3.0) {
425 426
		printf (
			"Error reading field origin: Read values (%f,%f,%f)\n",
gsell's avatar
gsell committed
427
			x_origin, y_origin, z_origin);
428
	}
gsell's avatar
gsell committed
429
	H5Block3dGetFieldSpacing (
430 431 432 433
		f, "TestField",
		&x_spacing,
		&y_spacing,
		&z_spacing );
gsell's avatar
gsell committed
434
	if (x_spacing != 2.0 || y_spacing != 3.0 || z_spacing != 4.0) {
435 436
		printf (
			"Error reading field spacing: Read values (%f,%f,%f)\n",
gsell's avatar
gsell committed
437
			x_spacing, y_spacing, z_spacing);
438 439
	}

gsell's avatar
gsell committed
440
	return H5_SUCCESS;
441 442 443 444 445 446 447 448 449 450 451 452
}

static h5_int64_t
_read_file (
	const char *fname,
	const int myproc,
	MPI_Comm comm,
	struct H5BlockPartition *layout
	) {
	printf ("PROC[%d]: Open file \"%s\" for reading ...\n",
		myproc, fname );

gsell's avatar
gsell committed
453 454 455 456 457 458 459
	h5_file_t file = H5OpenFile (fname, H5_O_RDONLY, H5_PROP_DEFAULT);
        H5SetStep (file, 0);

	_read_data (file, myproc, layout);
        _read_attributes (file, myproc, comm);

	H5CloseFile (file);
460
	
gsell's avatar
gsell committed
461
	return H5_SUCCESS;
462 463 464 465 466 467 468
}

int
main (
	int argc,
	char **argv
	) {
gsell's avatar
gsell committed
469
        h5_int64_t verbosity = DEFAULT_VERBOSITY;
470 471 472 473 474 475 476 477 478 479 480 481
	char *fname;
	int opt_with_ghosts = 0;
	int opt_read = 0;
	int opt_write = 0;
	struct H5BlockPartition *layout;

        if (argc == 1) {
                fprintf ( stderr,
                          "Usage: %s -w|-r [-g]\n",
                          argv[0] );
                return 1;
        }
gsell's avatar
gsell committed
482 483
	while (--argc) {
		if (strcmp (argv[argc], "-r") == 0)
484
			opt_read = 1;
gsell's avatar
gsell committed
485
		else if (strcmp (argv[argc], "-w") == 0)
486
			opt_write = 1;
gsell's avatar
gsell committed
487
		else if (strcmp (argv[argc], "-g") == 0)
488 489
			opt_with_ghosts = 1;
		else {
gsell's avatar
gsell committed
490 491 492 493
			fprintf (stderr,
				 "Illegal option %s\n\n"
				 "Usage: %s -w -r -g\n",
				 argv[argc], argv[0]);
494 495 496 497
			return 1;
		}
	}

gsell's avatar
gsell committed
498 499 500 501 502 503 504 505 506 507 508 509
        // initialize MPI & H5hut
        int comm_rank = 0;
        int comm_size = 1;
        MPI_Init (&argc, &argv);
        MPI_Comm comm = MPI_COMM_WORLD;
        MPI_Comm_rank (comm, &comm_rank);
        MPI_Comm_size (comm, &comm_size);

        H5AbortOnError ();
        H5SetVerbosityLevel (verbosity);

	switch (comm_size) {
510 511
	case 1:
		fname  = "blockfile1.h5";
gsell's avatar
gsell committed
512
		layout = &Layout1[comm_rank];
513 514
		break;
	case 8:
gsell's avatar
gsell committed
515
		if (opt_with_ghosts) {
516
			fname  = "blockfile8G.h5";
gsell's avatar
gsell committed
517
			layout = &Layout8G[comm_rank];
518 519
		} else {
			fname  = "blockfile8.h5";
gsell's avatar
gsell committed
520
			layout = &Layout8[comm_rank];
521 522 523
		}
		break;
	case 16:
gsell's avatar
gsell committed
524
		if (opt_with_ghosts) {
525
			fname  = "blockfile16G.h5";
gsell's avatar
gsell committed
526
			layout = &Layout16G[comm_rank];
527 528
		} else {
			fname  = "blockfile16.h5";
gsell's avatar
gsell committed
529
			layout = &Layout16[comm_rank];
530 531 532
		}
		break;
	case 32:
gsell's avatar
gsell committed
533
		if (opt_with_ghosts) {
534
			fname  = "blockfile32G.h5";
gsell's avatar
gsell committed
535
			layout = &Layout32G[comm_rank];
536 537
		} else {
			fname  = "blockfile32.h5";
gsell's avatar
gsell committed
538
			layout = &Layout32[comm_rank];
539 540 541
		}
		break;
	default:
542 543
		printf ( "Run this test on %d, %d, %d or %d processor(s)!\n",
			 1, 8, 16, 32);
544 545 546
		return 1;
	}

gsell's avatar
gsell committed
547 548 549 550
	if (opt_write) {
		_write_file (fname, comm_rank, comm, layout);
	} else if (opt_read) {
		_read_file (fname, comm_rank, comm, layout);
551 552 553 554 555
	}

	MPI_Finalize();
	return 0;
}