sandbox/acastillo/output_fields/xdmf/output_xdmf.h

    Light data is stored using eXtensible Markup Language (XML) to describe the data model and the data format. Heavy data is composed of large datasets stored using the Hierarchical Data Format HDF5. Additional features include support for a large number of objects, file compression, a parallel I/O implementation through the MPI-IO or MPI POSIX drivers.

    To use HDF5 we need to “declare” the library in the qcc system

    echo “@include <hdf5.h>” > $BASILISK/ast/std/hdf5.h echo “@include <hdf5_hl.h>” > $BASILISK/ast/std/hdf5_hl.h

    and declare the datatypes at the end of $BASILISK/ast/defaults.h

    echo “typedef hid_t, hsize_t, herr_t, H5L_info_t;” >> $BASILISK/ast/defaults.h

    see https://groups.google.com/g/basilisk-fr/c/CM270hBSfWo

    #pragma autolink -lhdf5
    #include <hdf5.h>
    
    
    // Macro to simplify dealing with slices
    #define shortcut_slice(n,_alpha )				        \
      double alpha = (_alpha - n.x*x - n.y*y - n.z*z)/Delta;\
      if (fabs(alpha) > 0.87)                               \
        continue;                                           \
    
    
    #include "output_xdmf_helpers.h"
    
    void write_xdmf_light_data(scalar *slist, vector *vlist, char * file_name, char * subname, int num_cells=0, int num_points=0, int dim=dimension){
    
        #if defined(_OPENMP)
        int num_omp = omp_get_max_threads();
        omp_set_num_threads(1);
        #endif
    
        // Open a file for writing
    	char name[111];
    	sprintf(name, "%s.xmf", subname);
        FILE * fp = fopen(name, "w");
    
    	/* Write the header of the xdmf file. */
        write_xdmf_header(fp, file_name);
    	write_xdmf_topology(fp, dim, num_cells, num_points, t);
        write_xdmf_attributes(fp, num_cells, slist, vlist);
        write_xdmf_footer(fp);
    
        // Close the file
    	fflush(fp);
    	fclose (fp);
    
        #if defined(_OPENMP)
        omp_set_num_threads(num_omp);
        #endif
    }
    
    trace void output_xmf(scalar *slist, vector *vlist, char *subname)
    {
    	hid_t acc_tpl1;	   // File access template
    	hid_t file_id;	   // HDF5 file ID
    	hid_t group_id;	   // HDF5 group ID
    	hsize_t count[2];  // Hyperslab selection parameters
    	hsize_t offset[2]; // Offset for hyperslab
    
    	char name[111];					 // Buffer for file name construction
    	sprintf(name, "%s.h5", subname); // Construct the HDF5 file name
    
    	// Define a scalar mask for periodic conditions
    	scalar per_mask[];
    	foreach ()
    	{
    		per_mask[] = 1.;
    	}
    	
    
    	// Obtain the number of points and cells and get a marker to reconstruct the topology
    	int num_points = 0, num_cells = 0, num_points_loc = 0, num_cells_loc = 0;
    	count_points_and_cells(&num_points, &num_cells, &num_points_loc, &num_cells_loc, per_mask);
    
    	// Calculate offsets for parallel I/O
    	int offset_points[npe()], offset_cells[npe()];
    	calculate_offsets(offset_points, offset_cells, num_points_loc, num_cells_loc, offset);
    
    	// Initialize marker for topology reconstruction
    	vertex scalar marker[];
    	initialize_marker(marker, offset);
    
    	// Write the light data
    	if (pid() == 0)
    		write_xdmf_light_data(slist, vlist, name, subname, num_cells, num_points);
    
    	// Write the heavy data
    	// Setup file access template with parallel I/O access
    	acc_tpl1 = H5Pcreate(H5P_FILE_ACCESS);
    	H5Pset_fapl_mpio(acc_tpl1, MPI_COMM_WORLD, MPI_INFO_NULL);
    
    	// Create a new HDF5 file collectively
    	file_id = H5Fcreate(name, H5F_ACC_TRUNC, H5P_DEFAULT, acc_tpl1);
    
    	// Release file-access template
    	H5Pclose(acc_tpl1);
    
    	// Populate and write the topology dataset
    	long *topo_dset;
    	populate_topo_dset(&topo_dset, num_cells_loc, offset_cells, count, offset, per_mask, marker);	
    	create_chunked_dataset(file_id, count, offset, "/Topology", num_cells, num_cells_loc, pow(2, dimension), topo_dset, H5T_NATIVE_LONG, num_cells/npe());
    	free(topo_dset);
    
    	// Create group for mesh geometry data
    	group_id = H5Gcreate(file_id, "Geometry", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    
    	// Populate and write the points dataset
    	double *points_dset;
    	populate_points_dset(&points_dset, num_points_loc, offset_points, count, offset);
    	create_chunked_dataset(group_id, count, offset, "/Geometry/Points", num_points, num_points_loc, 3, points_dset, H5T_NATIVE_DOUBLE, num_cells/npe());
    	free(points_dset);
    	H5Gclose(group_id);
    
    	// Create group for cell data
    	group_id = H5Gcreate(file_id, "Cells", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    
    	// Allocate memory and write scalar datasets
    	double *scalar_dset = (double *)malloc(num_cells_loc * sizeof(double));
    	for (scalar s in slist)
    	{
    		char substamp[1024];
    		sprintf(substamp, "/Cells/%s", s.name);
    		populate_scalar_dset(s, scalar_dset, num_cells_loc, offset_cells, count, offset, per_mask);
    		create_chunked_dataset(group_id, count, offset, substamp, num_cells, num_cells_loc, 1, scalar_dset, H5T_NATIVE_DOUBLE, num_cells/npe());
    	}
    	free(scalar_dset);
    
    	// Allocate memory and write vector datasets
    	double *vector_dset = (double *)malloc(num_cells_loc * 3 * sizeof(double));
    	for (vector v in vlist)
    	{
    		char substamp[1024];
    		sprintf(substamp, "/Cells/%s", v.x.name);
    		populate_vector_dset(v, vector_dset, num_cells_loc, offset_cells, count, offset, per_mask);
    		create_chunked_dataset(group_id, count, offset, substamp, num_cells, num_cells_loc, 3, vector_dset, H5T_NATIVE_DOUBLE, num_cells/npe());
    	}
    	free(vector_dset);
    	H5Gclose(group_id);
    
    	// Close HDF5 resources
    	H5Fclose(file_id);
    }
    
    
    #if dimension == 3
    trace void output_xmf_slice(scalar *slist, vector *vlist, char *subname, coord n = {0, 0, 1}, double _alpha = 0)
    {
    	hid_t acc_tpl1;	   // File access template
    	hid_t file_id;	   // HDF5 file ID
    	hid_t group_id;	   // HDF5 group ID
    	hsize_t count[2];  // Hyperslab selection parameters
    	hsize_t offset[2]; // Offset for hyperslab
    
    	char name[111];					 // Buffer for file name construction
    	sprintf(name, "%s.h5", subname); // Construct the HDF5 file name
    
    	// Define a scalar mask for periodic conditions
    	scalar per_mask[];
    	foreach () {
    		per_mask[] = 0.;
    		shortcut_slice(n, _alpha);
    		if (alpha > 0.){
    			#if EMBED
    			per_mask[] = cs[];
    			#else
    			per_mask[] = 1.;
    			#endif
    		}
    	}
    	
    
    	// Obtain the number of points and cells and get a marker to reconstruct the topology
    	int num_points = 0, num_cells = 0, num_points_loc = 0, num_cells_loc = 0;
    	count_points_and_cells_slice(&num_points, &num_cells, &num_points_loc, &num_cells_loc, per_mask, n, _alpha);
    
    	// Calculate offsets for parallel I/O
    	int offset_points[npe()], offset_cells[npe()];
    	calculate_offsets(offset_points, offset_cells, num_points_loc, num_cells_loc, offset);
    
    	// Initialize marker for topology reconstruction
    	vertex scalar marker[];
    	initialize_marker_slice(marker, offset, n, _alpha);
    
    	// Write the light data
    	if (pid() == 0)
    		write_xdmf_light_data(slist, vlist, name, subname, num_cells, num_points, dim=dimension-1);
    
    	// Write the heavy data
    	// Setup file access template with parallel I/O access
    	acc_tpl1 = H5Pcreate(H5P_FILE_ACCESS);
    	H5Pset_fapl_mpio(acc_tpl1, MPI_COMM_WORLD, MPI_INFO_NULL);
    
    	// Create a new HDF5 file collectively
    	file_id = H5Fcreate(name, H5F_ACC_TRUNC, H5P_DEFAULT, acc_tpl1);
    
    	// Release file-access template
    	H5Pclose(acc_tpl1);
    
    	// Populate and write the topology dataset
    	long *topo_dset;
    	populate_topo_dset_slice(&topo_dset, num_cells_loc, offset_cells, count, offset, per_mask, marker, n, _alpha);
    	create_chunked_dataset(file_id, count, offset, "/Topology", num_cells, num_cells_loc, pow(2, dimension-1), topo_dset, H5T_NATIVE_LONG, num_cells/npe());
    	free(topo_dset);
    
    	// Create group for mesh geometry data
    	group_id = H5Gcreate(file_id, "Geometry", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    
    	// Populate and write the points dataset
    	double *points_dset;
    	populate_points_dset_slice(&points_dset, num_points_loc, offset_points, count, offset, n, _alpha);
    	create_chunked_dataset(group_id, count, offset, "/Geometry/Points", num_points, num_points_loc, 3, points_dset, H5T_NATIVE_DOUBLE, num_points/npe());
    	free(points_dset);
    	H5Gclose(group_id);
    
    	// Create group for cell data
    	group_id = H5Gcreate(file_id, "Cells", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    
    	// Allocate memory and write scalar datasets
    	double *scalar_dset = (double *)malloc(num_cells_loc * sizeof(double));
    	for (scalar s in slist){
    		char substamp[1024];
    		sprintf(substamp, "/Cells/%s", s.name);
    		populate_scalar_dset_slice(s, scalar_dset, num_cells_loc, offset_cells, count, offset, per_mask, n, _alpha);
    		create_chunked_dataset(group_id, count, offset, substamp, num_cells, num_cells_loc, 1, scalar_dset, H5T_NATIVE_DOUBLE, num_cells/npe());
    	}
    	free(scalar_dset);
    
    	// Allocate memory and write vector datasets
    	double *vector_dset = (double *)malloc(num_cells_loc * 3 * sizeof(double));
    	for (vector v in vlist)	{
    		char substamp[1024];
    		sprintf(substamp, "/Cells/%s", v.x.name);
    		populate_vector_dset_slice(v, vector_dset, num_cells_loc, offset_cells, count, offset, per_mask, n, _alpha);
    		create_chunked_dataset(group_id, count, offset, substamp, num_cells, num_cells_loc, 3, vector_dset, H5T_NATIVE_DOUBLE, num_cells/npe());
    	}
    	free(vector_dset);
    	H5Gclose(group_id);
    
    	// Close HDF5 resources
    	H5Fclose(file_id);
    }
    #endif
    
    #undef shortcut_slice