Parallel I/O using Parallel HDF5

// written by Hoi Yen Loo, 1 September 2012 - 31 March 2013

#include <algorithm>

using std::string;
using std::stringstream;
using std::max;
using std::map;

#define TIME_LENGTH 10

// read in and write out structs for HDF treating
typedef struct {
	DT_INDEX init_grid_spac[3]; // initial grid spacings for top level grid
	DT_INDEX refine_grid_spac[3]; //  grid spacings for adaptive refinement
	DT_INDEX block_size[3]; // block size for each grid
	DT_VAL current_time;	// current time
	DT_VAL vis_time; // time interval in which to visualize data
	DT_VAL dt; // time discretisation step
	DT_VAL end_time; // end time of computation
	DT_VAL alpha; // heat diffusion coefficient
	DT_VAL mu; // dynamic viscosity of the fluid
	DT_VAL g[3]; // gravitational force
	DT_VAL err_tol; // maxmimal error tolerance in Jacoby-Step
	DT_INDEX max_iter; // maximal amount of iterations in Jacoby-Step
	DT_VAL omega; // SOR relaxation
	DT_TYPE write_vtk_result_array[CE_ARRAY_SIZE]; // vtk visualization selection
} fdomain_data;

typedef struct {
	DT_TAG mytag; // grid mytag
	DT_UID supergrid_uid; // uid of supergrid
	DT_DEPTH depth; // depth of the grid
	DT_TYPE ref_type; // indicate whether the grid is refined
	DT_INDEX n[3]; // block size of the grid
	DT_INDEX size_subgrids[3]; // size of the subgrid
} fgrid_data;

FHDF::FHDF(const MPI_Comm &_comm, const int &_comp_size, const int &_rank, FDomain *_domain, bool &success) :
		comm(_comm), comp_size(_comp_size), rank(_rank), domain(_domain), amount_nxyz_incl_ghost(0), first_write(true), latest_time(0) {

	#ifdef USE_HDF5
		if (H5open() < 0) {
			if (rank == 0) printf("Error-HDF5: The HDF5 library at rank %i could not be initialized!\n", rank);
			success = false;
			return;

		} else {
			// grid compound datatype
			hsize_t vector_dim[] = {3};
			hid_t vector_DT_INDEX_tid = H5Tarray_create(DT_INDEX_HDF, 1, vector_dim);

			grid_tid = H5Tcreate (H5T_COMPOUND, sizeof(FGrid));
			H5Tinsert(grid_tid, "mytag", HOFFSET(FGrid, mytag), DT_TAG_HDF);
			H5Tinsert(grid_tid, "supergrid_uid", HOFFSET(FGrid, supergrid_uid), DT_UID_HDF);
			H5Tinsert(grid_tid, "depth", HOFFSET(FGrid, depth), DT_DEPTH_HDF);
			H5Tinsert(grid_tid, "ref_type", HOFFSET(FGrid, ref_type), DT_TYPE_HDF);
			H5Tinsert(grid_tid, "n", HOFFSET(FGrid, n), vector_DT_INDEX_tid);
			H5Tinsert(grid_tid, "size_subgrids", HOFFSET(FGrid, size_subgrids), vector_DT_INDEX_tid);
			H5Tclose(vector_DT_INDEX_tid);

			domain->fhdf = this;
		}

	#endif

	success = true;
}

FHDF::~FHDF() {

	#ifdef USE_HDF5

		H5Tclose(grid_tid);
		H5Pclose(parallel_dataset_pid);

	#endif
}

void FHDF::setFilename(std::string _filename) {
	filename = _filename;
}

bool FHDF::writeHDF5(int *offset_array) {

	if (!fabs(domain->current_time - latest_time) < 1e-5) {
		if (!writeGridData(offset_array)) return false;

	} else if (first_write) {
		if (!writeDomainData(offset_array)) return false;
	}

	return true;

}

inline void FHDF::defineConstant() {

	#ifdef USE_HDF5

		amount_nxyz_incl_ghost = (domain->block_size[X] + 2) * (domain->block_size[Y] + 2) * (domain->block_size[Z] + 2);

		max_subgrid = max(domain->init_grid_spac[X] * domain->init_grid_spac[Y] * domain->init_grid_spac[Z],
						  domain->refine_grid_spac[X] * domain->refine_grid_spac[Y] * domain->refine_grid_spac[Z]);

	#endif

}

bool FHDF::createFile() {

	#ifdef USE_HDF5

		hid_t parallel_file_pid = H5Pcreate(H5P_FILE_ACCESS);
		H5Pset_fapl_mpio(parallel_file_pid, comm, MPI_INFO_NULL);
		file_fid = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, parallel_file_pid);
		H5Pclose(parallel_file_pid);

		parallel_dataset_pid = H5Pcreate(H5P_DATASET_XFER);
		H5Pset_dxpl_mpio(parallel_dataset_pid, H5FD_MPIO_COLLECTIVE);

		// write HDF5_FILE_VERSION and number of computation processes
		common_gid = H5Gcreate(file_fid, "/common", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);

		hid_t file_version_sid = H5Screate(H5S_SCALAR);
		hid_t file_version_aid = H5Acreate(file_fid, "HDF5_DATA_VERSION", DT_TYPE_HDF, file_version_sid, H5P_DEFAULT, H5P_DEFAULT);
		DT_TYPE file_version = HDF5_FILE_VERSION;
		if (H5Awrite(file_version_aid, DT_TYPE_HDF, &file_version) < 0) {
			if (rank == 0) printf("ERROR: %s: Could not write HDF5_DATA_VERSION to HDF5!\n", filename.c_str());
			return false;
		}
		H5Sclose(file_version_sid);
		H5Aclose(file_version_aid);

		hid_t comp_size_sid = H5Screate(H5S_SCALAR);
		hid_t comp_size_aid = H5Acreate(file_fid, "SIZE", DT_RANK_HDF, comp_size_sid, H5P_DEFAULT, H5P_DEFAULT);
		if (H5Awrite(comp_size_aid, DT_RANK_HDF, &comp_size) < 0) {
			if (rank == 0)printf("ERROR: %s: Could not write SIZE to HDF5!\n", filename.c_str());
			return false;
		}
		H5Sclose(comp_size_sid);
		H5Aclose(comp_size_aid);

	#endif

	return true;

}

bool FHDF::writeDomainData(const int *offset_array) {

	#ifndef USE_HDF5
		if (rank == 0) {
			printf("WARNING-WRITING HDF5: %s: No writing to %s because HDF5 capabilities is not activated!\n", filename.c_str(), filename.c_str() );
		}
		return true;
	#else

		createFile();

		hid_t domain_sid = H5Screate(H5S_SCALAR);

		// domain datatype
		hsize_t vector_dim[] = { 3 };
		hid_t domain_vector_DT_INDEX_tid = H5Tarray_create(DT_INDEX_HDF, 1, vector_dim);
		hid_t domain_array_DT_VAL_tid = H5Tarray_create(DT_VAL_HDF, 1, vector_dim);

		hsize_t array_vis_dim[] = { CE_ARRAY_SIZE };
		hid_t domain_array_DT_TYPE_tid = H5Tarray_create(DT_TYPE_HDF, 1, array_vis_dim);

		hid_t domain_tid = H5Tcreate (H5T_COMPOUND, sizeof(FDomain));
		H5Tinsert(domain_tid, "init_grid_spac", HOFFSET(FDomain, init_grid_spac), domain_vector_DT_INDEX_tid);
		H5Tinsert(domain_tid, "refine_grid_spac", HOFFSET(FDomain, refine_grid_spac), domain_vector_DT_INDEX_tid);
		H5Tinsert(domain_tid, "block_size", HOFFSET(FDomain, block_size), domain_vector_DT_INDEX_tid);
		H5Tinsert(domain_tid, "current_time", HOFFSET(FDomain, current_time), DT_VAL_HDF);
		H5Tinsert(domain_tid, "vis_time", HOFFSET(FDomain, vis_time), DT_VAL_HDF);
		H5Tinsert(domain_tid, "dt", HOFFSET(FDomain, dt), DT_VAL_HDF);
		H5Tinsert(domain_tid, "end_time", HOFFSET(FDomain, end_time), DT_VAL_HDF);
		H5Tinsert(domain_tid, "alpha", HOFFSET(FDomain, alpha), DT_VAL_HDF);
		H5Tinsert(domain_tid, "mu", HOFFSET(FDomain, mu), DT_VAL_HDF);
		H5Tinsert(domain_tid, "g", HOFFSET(FDomain, g), domain_array_DT_VAL_tid);
		H5Tinsert(domain_tid, "err_tol", HOFFSET(FDomain, err_tol), DT_VAL_HDF);
		H5Tinsert(domain_tid, "max_iter", HOFFSET(FDomain, max_iter), DT_INDEX_HDF);
		H5Tinsert(domain_tid, "omega", HOFFSET(FDomain, omega), DT_VAL_HDF);
		H5Tinsert(domain_tid, "write vtk result array", HOFFSET(FDomain, write_vtk_result_array), domain_array_DT_TYPE_tid);
		H5Tclose(domain_vector_DT_INDEX_tid);
		H5Tclose(domain_array_DT_VAL_tid);
		H5Tclose(domain_array_DT_TYPE_tid);

		hid_t domain_did = H5Dcreate(common_gid, "domain property", domain_tid, domain_sid, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);

		if (H5Dwrite(domain_did, domain_tid, domain_sid, H5S_ALL, H5P_DEFAULT, domain) < 0) {
			if (rank == 0) printf("ERROR: %s: Could not write domain data to HDF5!\n", filename.c_str());
			return false;
		}

		H5Tclose(domain_tid);
		H5Dclose(domain_did);
		H5Sclose(domain_sid);

		H5Gclose(common_gid);

		// create simulation group
		simulation_gid = H5Gcreate(file_fid, "/simulation", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);

		defineConstant();

		writeGridData(offset_array);

		return true;

	#endif

}

bool FHDF::writeGridData(const int *offset_array) {

	#ifdef USE_HDF5

		if (!first_write) {
			hid_t parallel_file_pid = H5Pcreate(H5P_FILE_ACCESS);
			H5Pset_fapl_mpio(parallel_file_pid, comm, MPI_INFO_NULL);
			file_fid = H5Fopen(filename.c_str(), H5F_ACC_RDWR, parallel_file_pid);
			H5Pclose(parallel_file_pid);

			simulation_gid = H5Gopen(file_fid, "/simulation", H5P_DEFAULT);
		}

		latest_time = domain->current_time;

		// write current time step for restarting simulation
		stringstream ss;
		ss << "/simulation/" << domain->current_time;
		string time_str = ss.str();

		// write latest time for restarting simulation
		ss.str("");
		ss << domain->current_time;

		char time_char[TIME_LENGTH];
		strcpy((char *)&time_char, ss.str().c_str());
		hid_t time_tid = H5Tcopy(H5T_C_S1);
		H5Tset_size(time_tid, TIME_LENGTH);

		hid_t time_sid;
		hid_t time_aid;
		if (first_write) {
			time_sid = H5Screate(H5S_SCALAR);
			time_aid = H5Acreate(simulation_gid, "TIME", time_tid, time_sid, H5P_DEFAULT, H5P_DEFAULT);
			first_write = false;

		} else {
			time_aid = H5Aopen(simulation_gid, "TIME", H5P_DEFAULT);
			time_sid = H5Aget_space(time_aid);
		}

		if (H5Awrite(time_aid, time_tid, time_char) < 0) {
			if (rank == 0) printf("ERROR: %s: Could not write TIME to HDF5!\n", filename.c_str());
			return false;
		}
		H5Tclose(time_tid);
		H5Sclose(time_sid);
		H5Aclose(time_aid);

		// write time group
		hid_t time_gid = H5Gcreate(file_fid, time_str.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);

		// write offset array for restarting simulation
		hsize_t offset_array_dim = comp_size + 2;
		hid_t offset_array_sid = H5Screate_simple(1, &offset_array_dim, NULL);
		hid_t offset_array_aid = H5Acreate(time_gid, "OFFSET ARRAY", DT_RANK_HDF, offset_array_sid, H5P_DEFAULT, H5P_DEFAULT);
		if (H5Awrite(offset_array_aid, DT_RANK_HDF, &(offset_array)[0]) < 0) {
			if (rank == 0) printf("ERROR: %s: Could not write OFFSET ARRAY to HDF5!\n", filename.c_str());
			return false;;
		}
		H5Sclose(offset_array_sid);
		H5Aclose(offset_array_aid);

		// write time in DT_VAL for setting domain->current_time
		hid_t time_DT_VAL_sid = H5Screate(H5S_SCALAR);
		hid_t time_DT_VAL_aid = H5Acreate(time_gid, "TIME", DT_VAL_HDF, time_DT_VAL_sid, H5P_DEFAULT, H5P_DEFAULT);
		H5Awrite(time_DT_VAL_aid, DT_VAL_HDF, &(domain->current_time));
		H5Sclose(time_DT_VAL_sid);
		H5Aclose(time_DT_VAL_aid);

		// write grid information
		hsize_t grid_dim[] = { offset_array[comp_size], 1 };
		hsize_t bbox_dim[] = { offset_array[comp_size], 6 };
		hsize_t subgrid_dim[] = { offset_array[comp_size], max_subgrid };
		hsize_t cell_data_dim[] = { offset_array[comp_size], amount_nxyz_incl_ghost * CE_ARRAY_SIZE };
		hsize_t temp_cell_data_dim[] = { offset_array[comp_size], amount_nxyz_incl_ghost * CT_TEMP_ARRAY_SIZE };
		hsize_t cell_type_dim[] = { offset_array[comp_size], amount_nxyz_incl_ghost };

		hid_t grid_sid = H5Screate_simple(2, grid_dim, NULL);
		hid_t bbox_sid = H5Screate_simple(2, bbox_dim, NULL);
		hid_t subgrid_sid = H5Screate_simple(2, subgrid_dim, NULL);
		hid_t current_cell_data_sid = H5Screate_simple(2, cell_data_dim, NULL);
		hid_t previous_cell_data_sid = H5Screate_simple(2, cell_data_dim, NULL);
		hid_t temp_cell_data_sid = H5Screate_simple(2, temp_cell_data_dim, NULL);
		hid_t cell_type_sid = H5Screate_simple(2, cell_type_dim, NULL);

		hid_t grid_did = H5Dcreate(time_gid, "grid property", grid_tid, grid_sid, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
		hid_t bbox_did = H5Dcreate(time_gid, "bounding box", DT_GEOM_HDF, bbox_sid, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
		hid_t subgrid_did = H5Dcreate(time_gid, "subgrid uid", DT_UID_HDF, subgrid_sid, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
		hid_t current_cell_data_did = H5Dcreate(time_gid, "current cell data", DT_VAL_HDF, current_cell_data_sid, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
		hid_t previous_cell_data_did = H5Dcreate(time_gid, "previous cell data", DT_VAL_HDF, previous_cell_data_sid, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
		hid_t temp_cell_data_did = H5Dcreate(time_gid, "temp cell data", DT_VAL_HDF, temp_cell_data_sid, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
		hid_t cell_type_did = H5Dcreate(time_gid, "cell type", DT_TYPE_HDF, cell_type_sid, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);

		grid_dim[0] = 1;
		bbox_dim[0] = 1;
		subgrid_dim[0] = 1;
		cell_data_dim[0] = 1;
		temp_cell_data_dim[0] = 1;
		cell_type_dim[0] = 1;
		hsize_t write_grid_offset[] = { offset_array[rank], 0 };

		int grid_each_rank = domain->all_grids.size();

		hid_t grid_memspace_sid;
		hid_t bbox_memspace_sid;
		hid_t subgrid_memspace_sid;
		hid_t current_cell_data_memspace_sid;
		hid_t previous_cell_data_memspace_sid;
		hid_t temp_cell_data_memspace_sid;
		hid_t cell_type_memspace_sid;

		hid_t grid_filespace_sid;
		hid_t bbox_filespace_sid;
		hid_t subgrid_filespace_sid;
		hid_t current_cell_data_filespace_sid;
		hid_t previous_cell_data_filespace_sid;
		hid_t temp_cell_data_filespace_sid;
		hid_t cell_type_filespace_sid;

		// -----------------------------------------------------
		// write grid data
		//------------------------------------------------------
		map<DT_GRIDID, FGrid*>::const_iterator grid_it = domain->all_grids.begin();
		for (int num_grid = 0; num_grid != offset_array[comp_size + 1]; ++num_grid) {

			grid_filespace_sid = H5Dget_space(grid_did);
			bbox_filespace_sid = H5Dget_space(bbox_did);
			current_cell_data_filespace_sid = H5Dget_space(current_cell_data_did);
			previous_cell_data_filespace_sid = H5Dget_space(previous_cell_data_did);
			temp_cell_data_filespace_sid = H5Dget_space(temp_cell_data_did);
			cell_type_filespace_sid = H5Dget_space(cell_type_did);

			FGrid *grid = 0;
			DT_GEOM *bbox = 0;
			DT_VAL *current_cell_data = 0;
			DT_VAL *previous_cell_data = 0;
			DT_VAL *temp_cell_data = 0;
			DT_TYPE *cell_type = 0;

			if (num_grid < grid_each_rank) {

				grid = &(grid_it->second)[0];
				bbox = &(grid_it->second->bbox)[0];
				current_cell_data = &(grid_it->second->array_val_current)[0];
				previous_cell_data = &(grid_it->second->array_val_previous)[0];
				temp_cell_data = &(grid_it->second->array_val_temp)[0];
				cell_type = &(grid_it->second->array_type)[0];

			} else {
				grid_dim[0] = grid_dim[1] = 0;
				bbox_dim[0] = bbox_dim[1] = 0;
				subgrid_dim[0] = subgrid_dim[1] = 0;
				cell_data_dim[0] = cell_data_dim[1] = 0;
				temp_cell_data_dim[0] = temp_cell_data_dim[1] = 0;
				cell_type_dim[0] = cell_type_dim[1] = 0;
				write_grid_offset[0] = write_grid_offset[1] = 0;
			}

			grid_memspace_sid = H5Screate_simple(2, grid_dim, NULL);
			bbox_memspace_sid = H5Screate_simple(2, bbox_dim, NULL);
			current_cell_data_memspace_sid = H5Screate_simple(2, cell_data_dim, NULL);
			previous_cell_data_memspace_sid = H5Screate_simple(2, cell_data_dim, NULL);
			temp_cell_data_memspace_sid = H5Screate_simple(2, temp_cell_data_dim, NULL);
			cell_type_memspace_sid = H5Screate_simple(2, cell_type_dim, NULL);

			H5Sselect_hyperslab(grid_filespace_sid, H5S_SELECT_SET, write_grid_offset, NULL, grid_dim, NULL);
			H5Sselect_hyperslab(bbox_filespace_sid, H5S_SELECT_SET, write_grid_offset, NULL, bbox_dim, NULL);
			H5Sselect_hyperslab(current_cell_data_filespace_sid, H5S_SELECT_SET, write_grid_offset, NULL, cell_data_dim, NULL);
			H5Sselect_hyperslab(previous_cell_data_filespace_sid, H5S_SELECT_SET, write_grid_offset, NULL, cell_data_dim, NULL);
			H5Sselect_hyperslab(temp_cell_data_filespace_sid, H5S_SELECT_SET, write_grid_offset, NULL, temp_cell_data_dim, NULL);
			H5Sselect_hyperslab(cell_type_filespace_sid, H5S_SELECT_SET, write_grid_offset, NULL, cell_type_dim, NULL);

			if (H5Dwrite(grid_did, grid_tid, grid_memspace_sid, grid_filespace_sid, parallel_dataset_pid, &grid[0]) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write grid property to HDF5!\n", filename.c_str());
				return false;
			}
			if (H5Dwrite(bbox_did, DT_GEOM_HDF, bbox_memspace_sid, bbox_filespace_sid, parallel_dataset_pid, bbox) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write bounding box to HDF5!\n", filename.c_str());
				return false;
			}
			if (H5Dwrite(current_cell_data_did, DT_VAL_HDF, current_cell_data_memspace_sid, current_cell_data_filespace_sid, parallel_dataset_pid, &current_cell_data[0]) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write current cell data to HDF5!\n", filename.c_str());
				return false;
			}
			if (H5Dwrite(previous_cell_data_did, DT_VAL_HDF, previous_cell_data_memspace_sid, previous_cell_data_filespace_sid, parallel_dataset_pid, &previous_cell_data[0]) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write previous cell data to HDF5!\n", filename.c_str());
				return false;
			}
			if (H5Dwrite(temp_cell_data_did, DT_VAL_HDF, temp_cell_data_memspace_sid, temp_cell_data_filespace_sid, parallel_dataset_pid, &temp_cell_data[0]) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write temp cell data to HDF5!\n", filename.c_str());
				return false;
			}
			if (H5Dwrite(cell_type_did, DT_TYPE_HDF, cell_type_memspace_sid, cell_type_filespace_sid, parallel_dataset_pid, &cell_type[0]) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write cell type to HDF5!\n", filename.c_str());
				return false;
			}

			H5Sclose(grid_filespace_sid);
			H5Sclose(grid_memspace_sid);
			H5Sclose(bbox_filespace_sid);
			H5Sclose(bbox_memspace_sid);
			H5Sclose(current_cell_data_filespace_sid);
			H5Sclose(current_cell_data_memspace_sid);
			H5Sclose(previous_cell_data_filespace_sid);
			H5Sclose(previous_cell_data_memspace_sid);
			H5Sclose(temp_cell_data_filespace_sid);
			H5Sclose(temp_cell_data_memspace_sid);
			H5Sclose(cell_type_filespace_sid);
			H5Sclose(cell_type_memspace_sid);

			// -----------------------------------------------------
			// write subgrid_uid
			//------------------------------------------------------
			DT_UID *write_subgrid = 0;

			int counter = 0;
			int subgrid_each_grid = 0;
			if (num_grid < grid_each_rank) {
				subgrid_each_grid = grid_it->second->subgrids_uid.size();
			}

			if(subgrid_each_grid != 0) {
				write_subgrid = (DT_UID *) calloc (max_subgrid, sizeof(DT_UID));

				DT_GRIDHASH ghash = 0;
				std::map<DT_GRIDHASH,DT_UID>::const_iterator sgit;
				for(DT_INDEX i = 0; i < grid_it->second->size_subgrids[X]; ++i) {
					for(DT_INDEX j = 0; j < grid_it->second->size_subgrids[Y]; ++j) {
						for(DT_INDEX k = 0; k < grid_it->second->size_subgrids[Z]; ++k) {
							ghash = UID::computeGridHash(i,j,k);

							sgit = grid_it->second->subgrids_uid.find(ghash);
							if(sgit == grid_it->second->subgrids_uid.end()) {
								printf("WARNING: Unable to find subgrid with hash %u\n", ghash);
							}

							write_subgrid[counter] = sgit->second;
							++counter;

						}
					}
				}
			}

			subgrid_filespace_sid = H5Dget_space(subgrid_did);

			if (subgrid_each_grid != 0) {
				subgrid_dim[0] = 1;
				subgrid_dim[1] = max_subgrid;

			} else {
				subgrid_dim[0] = subgrid_dim[1] = 0;
			}

			H5Sselect_hyperslab(subgrid_filespace_sid, H5S_SELECT_SET, write_grid_offset, NULL, subgrid_dim, NULL);

			subgrid_memspace_sid = H5Screate_simple(2, subgrid_dim, NULL);
			if (H5Dwrite(subgrid_did, DT_UID_HDF, subgrid_memspace_sid, subgrid_filespace_sid, parallel_dataset_pid, &(write_subgrid)[0]) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write subgrid uid to HDF5!\n", filename.c_str());
				return false;
			}

			H5Sclose(subgrid_filespace_sid);
			H5Sclose(subgrid_memspace_sid);

			if (write_subgrid != 0) free(write_subgrid);

			// -----------------------------------------------------
			// finish writing subgrid_uid for a particular grid
			// -----------------------------------------------------

			++write_grid_offset[0];
			++grid_it;

		}
		// -----------------------------------------------------
		// finish writing grid data
		//------------------------------------------------------

		H5Dclose(grid_did);
		H5Sclose(grid_sid);

		H5Dclose(bbox_did);
		H5Sclose(bbox_sid);

		H5Dclose(subgrid_did);
		H5Sclose(subgrid_sid);

		H5Dclose(current_cell_data_did);
		H5Sclose(current_cell_data_sid);

		H5Dclose(previous_cell_data_did);
		H5Sclose(previous_cell_data_sid);

		H5Dclose(temp_cell_data_did);
		H5Sclose(temp_cell_data_sid);

		H5Dclose(cell_type_did);
		H5Sclose(cell_type_sid);

		H5Gclose(time_gid);
		H5Gclose(simulation_gid);
		H5Fclose(file_fid);

	#endif

	return true;
}

bool FHDF::readHDF5() {

	#ifndef USE_HDF5

		if (rank == 0) {
			printf("ERROR-READING HDF5: %s: No reading from %s because HDF5 capabilities is not activated!\n",
					filename.c_str(), filename.c_str() );
		}
	return false;
	#else

		hid_t parallel_file_pid = H5Pcreate(H5P_FILE_ACCESS);
		H5Pset_fapl_mpio(parallel_file_pid, comm, MPI_INFO_NULL);
		hid_t read_file_fid = H5Fopen(filename.c_str(), H5F_ACC_RDWR, parallel_file_pid);
		H5Pclose(parallel_file_pid);

		parallel_dataset_pid = H5Pcreate(H5P_DATASET_XFER);
		H5Pset_dxpl_mpio(parallel_dataset_pid, H5FD_MPIO_COLLECTIVE);

		// check file version
		hid_t file_version_aid = H5Aopen(read_file_fid, "HDF5_DATA_VERSION", H5P_DEFAULT);
		DT_TYPE file_version;
		H5Aread(file_version_aid, DT_TYPE_HDF, &file_version);
		H5Aclose(file_version_aid);

		if (file_version != HDF5_FILE_VERSION) {
			if (rank == 0) {
				printf("ERROR-READING HDF5: %s: Mismatch between file version of %s and HDF5_FILE_VERSION of the code!\n"
						"      Please rerun the simulation with the correct file version!\n", filename.c_str(), filename.c_str() );
			}
			return false;
		}

		// check number of processes
		hid_t size_aid = H5Aopen(read_file_fid, "SIZE", H5P_DEFAULT);
		DT_RANK size;
		H5Aread(size_aid, DT_RANK_HDF, &size);
		H5Aclose(size_aid);

		if (size != comp_size) {
			if (rank == 0) {
				printf("ERROR-READING HDF5: %s: Mismatch between number of available processes and number of processes required for the grid distribution!\n"
						"      Please rerun the simulation with %i computation processes!\n", filename.c_str(), size );
			}
			return false;
		}

		// -----------------------------------------------------
		// read domain data
		// -----------------------------------------------------

		hsize_t vector_dim[] = { 3 };
		hid_t domain_vector_DT_INDEX_tid = H5Tarray_create(DT_INDEX_HDF, 1, vector_dim);
		hid_t domain_array_DT_VAL_tid = H5Tarray_create(DT_VAL_HDF, 1, vector_dim);

		hsize_t array_vis_dim[] = { CE_ARRAY_SIZE };
		hid_t domain_array_DT_TYPE_tid = H5Tarray_create(DT_TYPE_HDF, 1, array_vis_dim);

		hid_t domain_tid = H5Tcreate (H5T_COMPOUND, sizeof(fdomain_data));
		H5Tinsert(domain_tid, "init_grid_spac", HOFFSET(fdomain_data, init_grid_spac), domain_vector_DT_INDEX_tid);
		H5Tinsert(domain_tid, "refine_grid_spac", HOFFSET(fdomain_data, refine_grid_spac), domain_vector_DT_INDEX_tid);
		H5Tinsert(domain_tid, "block_size", HOFFSET(fdomain_data, block_size), domain_vector_DT_INDEX_tid);
		H5Tinsert(domain_tid, "current_time", HOFFSET(fdomain_data, current_time), DT_VAL_HDF);
		H5Tinsert(domain_tid, "vis_time", HOFFSET(fdomain_data, vis_time), DT_VAL_HDF);
		H5Tinsert(domain_tid, "dt", HOFFSET(fdomain_data, dt), DT_VAL_HDF);
		H5Tinsert(domain_tid, "end_time", HOFFSET(fdomain_data, end_time), DT_VAL_HDF);
		H5Tinsert(domain_tid, "alpha", HOFFSET(fdomain_data, alpha), DT_VAL_HDF);
		H5Tinsert(domain_tid, "mu", HOFFSET(fdomain_data, mu), DT_VAL_HDF);
		H5Tinsert(domain_tid, "g", HOFFSET(fdomain_data, g), domain_array_DT_VAL_tid);
		H5Tinsert(domain_tid, "err_tol", HOFFSET(fdomain_data, err_tol), DT_VAL_HDF);
		H5Tinsert(domain_tid, "max_iter", HOFFSET(fdomain_data, max_iter), DT_INDEX_HDF);
		H5Tinsert(domain_tid, "omega", HOFFSET(fdomain_data, omega), DT_VAL_HDF);
		H5Tinsert(domain_tid, "write vtk result array", HOFFSET(fdomain_data, write_vtk_result_array), domain_array_DT_TYPE_tid);
		H5Tclose(domain_vector_DT_INDEX_tid);
		H5Tclose(domain_array_DT_VAL_tid);
		H5Tclose(domain_array_DT_TYPE_tid);

		hid_t read_common_gid = H5Gopen(read_file_fid, "/common", H5P_DEFAULT);
		hid_t read_domain_did = H5Dopen(read_common_gid, "domain property", H5P_DEFAULT);
		hid_t read_domain_sid = H5Dget_space(read_domain_did);

		fdomain_data *read_domain = new fdomain_data;
		H5Dread(read_domain_did, domain_tid, H5S_ALL, H5S_ALL, H5P_DEFAULT, read_domain);

		H5Dclose(read_domain_did);
		H5Sclose(read_domain_sid);
		H5Gclose(read_common_gid);

		for (int i = 0; i != 3; ++i) {
			domain->init_grid_spac[i] = read_domain->init_grid_spac[i];
			domain->refine_grid_spac[i] = read_domain->refine_grid_spac[i];
			domain->block_size[i] = read_domain->block_size[i];
			domain->g[i] = read_domain->g[i];
		}

		for (int i = 0; i != CE_ARRAY_SIZE; ++i) {
			domain->write_vtk_result_array[i] = read_domain->write_vtk_result_array[i];
		}

		domain->vis_time = read_domain->vis_time;
		domain->dt = read_domain->dt;
		domain->end_time = read_domain->end_time;
		domain->alpha = read_domain->alpha;
		domain->mu = read_domain->mu;
		domain->err_tol = read_domain->err_tol;
		domain->max_iter = read_domain->max_iter;
		domain->omega = read_domain->omega;

		if (rank == 0) domain->domain_with_master_grid = true;
		else domain->domain_with_master_grid = false;

		// -----------------------------------------------------
		// read TIME
		// -----------------------------------------------------

		hid_t read_simulation_gid = H5Gopen(read_file_fid, "/simulation", H5P_DEFAULT);

		hid_t time_aid = H5Aopen(read_simulation_gid, "TIME", H5P_DEFAULT);

		char read_time[TIME_LENGTH];
		hid_t time_tid = H5Tcopy(H5T_C_S1);
		H5Tset_size(time_tid, TIME_LENGTH);

		H5Aread(time_aid, time_tid, read_time);
		H5Aclose(time_aid);
		H5Tclose(time_tid);

		// -----------------------------------------------------
		// read OFFSET ARRAY and DT_VAL TIME
		// -----------------------------------------------------

		string time = read_time;
		string time_str = "/simulation/" + time;

		hid_t read_time_gid = H5Gopen(read_file_fid, time_str.c_str(), H5P_DEFAULT);

		hid_t offset_array_aid = H5Aopen(read_time_gid, "OFFSET ARRAY", H5P_DEFAULT);
		DT_RANK *offset_array = new DT_RANK[comp_size + 2];
		H5Aread(offset_array_aid, DT_RANK_HDF, &(offset_array)[0]);
		H5Aclose(offset_array_aid);

		int grid_per_rank = offset_array[rank + 1] - offset_array[rank];

		hid_t time_DT_VAL_aid = H5Aopen(read_time_gid, "TIME", H5P_DEFAULT);
		H5Aread(time_DT_VAL_aid, DT_VAL_HDF, &(domain->current_time));
		H5Aclose(time_DT_VAL_aid);

		read_domain->current_time = domain->current_time;

		if (fabs(domain->current_time - domain->end_time) < 1e-8) {
			if (rank == 0) {
				printf("ERROR: %s: No simulation because the current time is equal to the end time!\n"
						"      Please change the end time in %s and rerun the simulation!\n", filename.c_str(), filename.c_str() );
			}
			return false;
		}

		// -----------------------------------------------------
		// write domain data
		// -----------------------------------------------------

		filename = filename.substr(0, filename.size() - 3) + "_restart.h5";
		createFile();

		hid_t domain_sid = H5Screate(H5S_SCALAR);
		hid_t domain_did = H5Dcreate(file_fid, "domain property", domain_tid, domain_sid, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);

		H5Dwrite(domain_did, domain_tid, domain_sid, H5S_ALL, H5P_DEFAULT, read_domain);

		H5Tclose(domain_tid);
		H5Sclose(domain_sid);
		H5Dclose(domain_did);
		H5Gclose(common_gid);

		delete read_domain;

		// -----------------------------------------------------
		// read grid data
		// -----------------------------------------------------

		hid_t read_grid_did = H5Dopen(read_time_gid, "grid property", H5P_DEFAULT);
		hid_t read_grid_filespace_sid = H5Dget_space(read_grid_did);
		hid_t read_bbox_did = H5Dopen(read_time_gid, "bounding box", H5P_DEFAULT);
		hid_t read_bbox_filespace_sid = H5Dget_space(read_bbox_did);
		hid_t read_subgrid_did = H5Dopen(read_time_gid, "subgrid uid", H5P_DEFAULT);
		hid_t read_subgrid_filespace_sid = H5Dget_space(read_subgrid_did);
		hid_t read_current_cell_data_did = H5Dopen(read_time_gid, "current cell data", H5P_DEFAULT);
		hid_t read_current_cell_data_filespace_sid = H5Dget_space(read_current_cell_data_did);
		hid_t read_previous_cell_data_did = H5Dopen(read_time_gid, "previous cell data", H5P_DEFAULT);
		hid_t read_previous_cell_data_filespace_sid = H5Dget_space(read_previous_cell_data_did);
		hid_t read_temp_cell_data_did = H5Dopen(read_time_gid, "temp cell data", H5P_DEFAULT);
		hid_t read_temp_cell_data_filespace_sid = H5Dget_space(read_temp_cell_data_did);
		hid_t read_cell_type_did = H5Dopen(read_time_gid, "cell type", H5P_DEFAULT);
		hid_t read_cell_type_filespace_sid = H5Dget_space(read_cell_type_did);

		hid_t vector_DT_INDEX_tid = H5Tarray_create(DT_INDEX_HDF, 1, vector_dim);

		hid_t read_grid_tid = H5Tcreate (H5T_COMPOUND, sizeof(fgrid_data));
		H5Tinsert(read_grid_tid, "mytag", HOFFSET(fgrid_data, mytag), DT_TAG_HDF);
		H5Tinsert(read_grid_tid, "supergrid_uid", HOFFSET(fgrid_data, supergrid_uid), DT_UID_HDF);
		H5Tinsert(read_grid_tid, "depth", HOFFSET(fgrid_data, depth), DT_DEPTH_HDF);
		H5Tinsert(read_grid_tid, "ref_type", HOFFSET(fgrid_data, ref_type), DT_TYPE_HDF);
		H5Tinsert(read_grid_tid, "n", HOFFSET(fgrid_data, n), vector_DT_INDEX_tid);
		H5Tinsert(read_grid_tid, "size_subgrids", HOFFSET(fgrid_data, size_subgrids), vector_DT_INDEX_tid);
		H5Tclose(vector_DT_INDEX_tid);

		defineConstant();

		hsize_t read_grid_dim[] = { 1, 1 };
		hsize_t read_bbox_dim[] = { 1, 6 };
		hsize_t read_subgrid_dim[] = { 1, max_subgrid };
		hsize_t read_cell_data_dim[] = { 1, amount_nxyz_incl_ghost * CE_ARRAY_SIZE };
		hsize_t read_temp_cell_data_dim[] = { 1, amount_nxyz_incl_ghost * CT_TEMP_ARRAY_SIZE };
		hsize_t read_cell_type_dim[] = { 1, amount_nxyz_incl_ghost };
		hsize_t read_grid_offset[] = { offset_array[rank], 0 };

		hid_t read_grid_memspace_sid = H5Screate_simple(2, read_grid_dim, NULL);
		hid_t read_bbox_memspace_sid = H5Screate_simple(2, read_bbox_dim, NULL);
		hid_t read_subgrid_memspace_sid = H5Screate_simple(2, read_subgrid_dim, NULL);
		hid_t read_current_cell_data_memspace_sid = H5Screate_simple(2, read_cell_data_dim, NULL);
		hid_t read_previous_cell_data_memspace_sid = H5Screate_simple(2, read_cell_data_dim, NULL);
		hid_t read_temp_cell_data_memspace_sid = H5Screate_simple(2, read_temp_cell_data_dim, NULL);
		hid_t read_cell_type_memspace_sid = H5Screate_simple(2, read_cell_type_dim, NULL);

		fgrid_data *read_grid = new fgrid_data;
		DT_GEOM *read_bbox = new DT_GEOM[6];
		DT_UID *read_subgrid = new DT_UID[max_subgrid];

		for (int num_grid = 0; num_grid != grid_per_rank; ++num_grid) {

			read_grid_dim[0] = read_grid_dim[1] = 1;
			read_grid_offset[0] = offset_array[rank] + num_grid;

			H5Sselect_hyperslab(read_grid_filespace_sid, H5S_SELECT_SET, read_grid_offset, NULL, read_grid_dim, NULL);
			H5Sselect_hyperslab(read_bbox_filespace_sid, H5S_SELECT_SET, read_grid_offset, NULL, read_bbox_dim, NULL);
			H5Sselect_hyperslab(read_subgrid_filespace_sid, H5S_SELECT_SET, read_grid_offset, NULL, read_subgrid_dim, NULL);
			H5Sselect_hyperslab(read_current_cell_data_filespace_sid, H5S_SELECT_SET, read_grid_offset, NULL, read_cell_data_dim, NULL);
			H5Sselect_hyperslab(read_previous_cell_data_filespace_sid, H5S_SELECT_SET, read_grid_offset, NULL, read_cell_data_dim, NULL);
			H5Sselect_hyperslab(read_temp_cell_data_filespace_sid, H5S_SELECT_SET, read_grid_offset, NULL, read_temp_cell_data_dim, NULL);
			H5Sselect_hyperslab(read_cell_type_filespace_sid, H5S_SELECT_SET, read_grid_offset, NULL, read_cell_type_dim, NULL);

			if (H5Dread(read_grid_did, read_grid_tid, read_grid_memspace_sid, read_grid_filespace_sid, H5P_DEFAULT, read_grid) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write grid property to HDF5!\n", filename.c_str());
				return false;
			}
			if (H5Dread(read_bbox_did, DT_GEOM_HDF, read_bbox_memspace_sid, read_bbox_filespace_sid, H5P_DEFAULT, &(read_bbox)[0]) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write bounding box to HDF5!\n", filename.c_str());
				return false;
			}
			if (H5Dread(read_subgrid_did, DT_UID_HDF, read_subgrid_memspace_sid, read_subgrid_filespace_sid, H5P_DEFAULT, &(read_subgrid)[0]) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write subgrid uid to HDF5!\n", filename.c_str());
				return false;
			}

			// construct new grid
			DT_GRIDID grid_id = UID::getGridIDFromTag(read_grid->mytag);
			DT_GRIDHASH grid_hash = UID::getGridHashFromTag(read_grid->mytag);
			FGrid* newgrid = new FGrid(domain, read_grid->supergrid_uid, grid_hash, read_bbox, read_grid->depth, grid_id);

			if(newgrid->n[X] != read_grid->n[X] || newgrid->n[Y] != read_grid->n[Y] || newgrid->n[Z] != read_grid->n[Z]) {
				printf("The fields n do not fit from rereading!!! Danger Danger!\n");
			}

			if(newgrid->size_subgrids[X] != read_grid->size_subgrids[X] || newgrid->size_subgrids[Y] != read_grid->size_subgrids[Y] || newgrid->size_subgrids[Z] != read_grid->size_subgrids[Z]) {
				printf("The fields n do not fit from rereading!!! Danger Danger!\n");
			}

			newgrid->subgrids_uid.clear();

			// set grid ref_type
			newgrid->ref_type = read_grid->ref_type;

			// construct subgrid
			int counter = 0;
			if(read_grid->ref_type == TP_GDREF_YES) {

				DT_GRIDHASH ghash = 0;

				for(DT_INDEX i = 0; i < newgrid->size_subgrids[X]; ++i) {
					for(DT_INDEX j = 0; j < newgrid->size_subgrids[Y]; ++j) {
						for(DT_INDEX k = 0; k < newgrid->size_subgrids[Z]; ++k) {
							ghash = UID::computeGridHash(i,j,k);

							newgrid->subgrids_uid[ghash] = read_subgrid[counter];
							++counter;
						}
					}
				}
			}

			// read current cell data, previous cell data, temp cell data and cell type
			if (H5Dread(read_current_cell_data_did, DT_VAL_HDF, read_current_cell_data_memspace_sid, read_current_cell_data_filespace_sid, H5P_DEFAULT, &(newgrid->array_val_current)[0]) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write current cell data to HDF5!\n", filename.c_str());
				return false;
			}
			if (H5Dread(read_previous_cell_data_did, DT_VAL_HDF, read_previous_cell_data_memspace_sid, read_previous_cell_data_filespace_sid, H5P_DEFAULT, &(newgrid->array_val_previous)[0]) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write previous cell data to HDF5!\n", filename.c_str());
				return false;
			}
			if (H5Dread(read_temp_cell_data_did, DT_VAL_HDF, read_temp_cell_data_memspace_sid, read_temp_cell_data_filespace_sid, H5P_DEFAULT, &(newgrid->array_val_temp)[0]) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write temp cell data to HDF5!\n", filename.c_str());
				return false;
			}
			if (H5Dread(read_cell_type_did, DT_TYPE_HDF, read_cell_type_memspace_sid, read_cell_type_filespace_sid, H5P_DEFAULT, &(newgrid->array_type)[0]) < 0) {
				if (rank == 0) printf("ERROR: %s: Could not write cell type to HDF5!\n", filename.c_str());
				return false;
			}

		}

		delete read_grid;
		delete [] read_bbox;
		delete [] read_subgrid;

		H5Tclose(read_grid_tid);
		H5Dclose(read_grid_did);
		H5Sclose(read_grid_filespace_sid);
		H5Sclose(read_grid_memspace_sid);

		H5Dclose(read_bbox_did);
		H5Sclose(read_bbox_filespace_sid);
		H5Sclose(read_bbox_memspace_sid);

		H5Dclose(read_subgrid_did);
		H5Sclose(read_subgrid_filespace_sid);
		H5Sclose(read_subgrid_memspace_sid);

		H5Dclose(read_current_cell_data_did);
		H5Sclose(read_current_cell_data_filespace_sid);
		H5Sclose(read_current_cell_data_memspace_sid);

		H5Dclose(read_previous_cell_data_did);
		H5Sclose(read_previous_cell_data_filespace_sid);
		H5Sclose(read_previous_cell_data_memspace_sid);

		H5Dclose(read_temp_cell_data_did);
		H5Sclose(read_temp_cell_data_filespace_sid);
		H5Sclose(read_temp_cell_data_memspace_sid);

		H5Dclose(read_cell_type_did);
		H5Sclose(read_cell_type_filespace_sid);
		H5Sclose(read_cell_type_memspace_sid);

		H5Gclose(read_time_gid);
		H5Gclose(read_simulation_gid);
		H5Fclose(read_file_fid);

		// -----------------------------------------------------
		// write grid data
		// -----------------------------------------------------

		simulation_gid = H5Gcreate(file_fid, "/simulation", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
		writeGridData(offset_array);

		return true;

	#endif

}