write_blocks_hdf5_r8.c

做网格的好程序

#include "hdf5.h"
#include "mpi.h"
#include <ctype.h>

#include "underscore.h"
#define float double
#define H5T_NATIVE_FLOAT H5T_NATIVE_DOUBLE

void CREATE_File_hdf5_r8 (char* file_name_in, hid_t* file_id) {

  hid_t plist_id;
  char file_name[80];
  int i; int ii;

  /*
     Don't know why I had to add this code, but it appears that some
     non-numeric, non-chararcter data gets inserted into this file name
     when called from the c_interface.  This code will remove any weird
     characters
  */
  for(i=0; i<80; i++) {
    if (!isalnum(file_name_in[i]) && !ispunct(file_name_in[i])) {
      for (ii=i;ii<79;ii++) {
        file_name_in[ii]=file_name_in[ii+1];
      }
    }
  }

  /* strip off fortran junk from the filename */
  strncpy(file_name, file_name_in, 80);
  for(i=79; i>=0; i--) {
    if (file_name[i] != ' ') { file_name[i+1]=0; break; }
  }

  /* Set up file access property list with parllel I/O access */
  plist_id = H5Pcreate(H5P_FILE_ACCESS);
  H5Pset_fapl_mpio(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL);

  /* Create a new file collectively */
  *file_id = H5Fcreate(file_name, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
  H5Pclose(plist_id);

}



void WRITE_Header_hdf5_r8(hid_t* file_id, int* tot_blocks,
			  float* user_attr_1,
			  float* user_attr_2,
			  float* user_attr_3,
			  float* user_attr_4,
			  float* user_attr_5) {

  hsize_t  dims_header[1];
  hid_t    dset_header, dataspace_header, group_id, attribute_id;
  herr_t   status;
  int      mpi_rank;
  int      temp[1];

  dims_header[0] = 1;
  dataspace_header = H5Screate_simple(1, dims_header, NULL);
  dset_header = H5Dcreate(*file_id, "Global Data", H5T_NATIVE_INT, 
			  dataspace_header, H5P_DEFAULT);

  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);

  if (mpi_rank == 0) {
    temp[0] = *tot_blocks;
    status = H5Dwrite(dset_header, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, 
                      H5P_DEFAULT, temp);
  }

  H5Sclose(dataspace_header);
  H5Dclose(dset_header);

  /* ADD USER ATTRIBUTES */
  group_id = H5Gopen(*file_id, "/");

  /* attribute 1 */
  dataspace_header = H5Screate(H5S_SCALAR);
  attribute_id = H5Acreate(group_id, "user attribute 1", H5T_NATIVE_FLOAT,
			   dataspace_header, H5P_DEFAULT);
  if (mpi_rank == 0) {
    status = H5Awrite(attribute_id, H5T_NATIVE_FLOAT, user_attr_1);
  }
  H5Sclose(dataspace_header);
  H5Aclose(attribute_id);

  /* attribute 2 */
  dataspace_header = H5Screate(H5S_SCALAR);
  attribute_id = H5Acreate(group_id, "user attribute 2", H5T_NATIVE_FLOAT,
			   dataspace_header, H5P_DEFAULT);
  if (mpi_rank == 0) {
    status = H5Awrite(attribute_id, H5T_NATIVE_FLOAT, user_attr_2);
  }
  H5Sclose(dataspace_header);
  H5Aclose(attribute_id);

  /* attribute 3 */
  dataspace_header = H5Screate(H5S_SCALAR);
  attribute_id = H5Acreate(group_id, "user attribute 3", H5T_NATIVE_FLOAT,
			   dataspace_header, H5P_DEFAULT);
  if (mpi_rank == 0) {
    status = H5Awrite(attribute_id, H5T_NATIVE_FLOAT, user_attr_3);
  }
  H5Sclose(dataspace_header);
  H5Aclose(attribute_id);

  /* attribute 4 */
  dataspace_header = H5Screate(H5S_SCALAR);
  attribute_id = H5Acreate(group_id, "user attribute 4", H5T_NATIVE_FLOAT,
			   dataspace_header, H5P_DEFAULT);
  if (mpi_rank == 0) {
    status = H5Awrite(attribute_id, H5T_NATIVE_FLOAT, user_attr_4);
  }
  H5Sclose(dataspace_header);
  H5Aclose(attribute_id);

  /* attribute 5 */
  dataspace_header = H5Screate(H5S_SCALAR);
  attribute_id = H5Acreate(group_id, "user attribute 5", H5T_NATIVE_FLOAT,
			   dataspace_header, H5P_DEFAULT);
  if (mpi_rank == 0) {
    status = H5Awrite(attribute_id, H5T_NATIVE_FLOAT, user_attr_5);
  }
  H5Sclose(dataspace_header);
  H5Aclose(attribute_id);

  H5Gclose(group_id);
      
}

 
void CREATE_Block_data_type_hdf5_r8(hid_t* block_type_id, 
				   int* mdim, 
				   int* ndim, 
				   int* ngid, 
				   int* mflags, 
				   int* nvar_chk_cc, 
				   int* nvar_chk_fc, 
				   int* nvar_chk_ec, 
				   int* nvar_chk_nc, 
				   int* il0, int* iu0, 
				   int* jl0, int* ju0, 
				   int* kl0, int* ku0) 
{
  
  int rank_array;
  hsize_t dims_array[1];
  int k2d, k3d;

  hid_t gid_tid, bflags_tid, coord_tid, bnd_box_tid, unk_tid;
  hid_t facevarx_tid, facevary_tid, facevarz_tid;
  hid_t unk_e_x_tid, unk_e_y_tid, unk_e_z_tid;
  hid_t unk_n_tid;
  int nx, ny, nz;
  int size_block_type;
  int offset;
  
  
  k2d = 0;
  if (*ndim >= 2) k2d = 1;
  k3d = 0;
  if (*ndim == 3) k3d = 1;

  /* Create the memory data type */

  size_block_type = (sizeof(float)*3) + 
    (sizeof(float) * *ngid) * + (sizeof(float) * *mflags) +
    (sizeof(float) * *mdim) + (sizeof(float) * 2 * *mdim) +
    (sizeof(float)) + 
    (sizeof(float) * *nvar_chk_cc * (*iu0-*il0) * (*ju0-*jl0) * (*ku0-*kl0)) +
    (sizeof(float) * *nvar_chk_fc * (*iu0-*il0+1) * (*ju0-*jl0) * (*ku0-*kl0)) +
    (sizeof(float) * *nvar_chk_fc * (*iu0-*il0) * (*ju0-*jl0+k2d) * (*ku0-*kl0)) +
    (sizeof(float) * *nvar_chk_fc * (*iu0-*il0) * (*ju0-*jl0) * (*ku0-*kl0+k3d)) +
    (sizeof(float) * *nvar_chk_ec * (*iu0-*il0) * (*ju0-*jl0+k2d) * (*ku0-*kl0+k3d)) +
    (sizeof(float) * *nvar_chk_ec * (*iu0-*il0+1) * (*ju0-*jl0) * (*ku0-*kl0+k3d)) +
    (sizeof(float) * *nvar_chk_ec * (*iu0-*il0+1) * (*ju0-*jl0+k2d) * (*ku0-*kl0)) +
    (sizeof(float) * *nvar_chk_nc * (*iu0-*il0+1) * (*ju0-*jl0+k2d) * (*ku0-*kl0+k3d));

  *block_type_id = H5Tcreate (H5T_COMPOUND, size_block_type);

  rank_array = 1;
  dims_array[0] = *ngid;
  gid_tid = H5Tarray_create(H5T_NATIVE_FLOAT, rank_array, 
			    dims_array, NULL);
  rank_array = 1;
  dims_array[0] = *mflags;
  bflags_tid = H5Tarray_create(H5T_NATIVE_FLOAT, rank_array, 
			       dims_array, NULL);
  rank_array = 1;
  dims_array[0] = *mdim;
  coord_tid = H5Tarray_create(H5T_NATIVE_FLOAT, rank_array, 
                              dims_array, NULL);
  rank_array = 1;
  dims_array[0] = 2 * *mdim;
  bnd_box_tid = H5Tarray_create(H5T_NATIVE_FLOAT, rank_array, 
				dims_array, NULL);
  rank_array = 1;
  nx = *iu0 - *il0;
  ny = *ju0 - *jl0;
  nz = *ku0 - *kl0;

  if (*nvar_chk_cc > 0) {
  dims_array[0] = *nvar_chk_cc * nx * ny * nz;
  unk_tid = H5Tarray_create(H5T_NATIVE_FLOAT, rank_array, 
			    dims_array, NULL);
  }

  if (*nvar_chk_fc > 0) {
  dims_array[0] = *nvar_chk_fc * (nx+1) * ny * nz;
  facevarx_tid = H5Tarray_create(H5T_NATIVE_FLOAT, rank_array, 
				 dims_array, NULL);
  dims_array[0] = *nvar_chk_fc * nx * (ny+k2d) * nz;
  facevary_tid = H5Tarray_create(H5T_NATIVE_FLOAT, rank_array, 
				 dims_array, NULL);
  dims_array[0] = *nvar_chk_fc * nx * ny * (nz+k3d);
  facevarz_tid = H5Tarray_create(H5T_NATIVE_FLOAT, rank_array, 
				 dims_array, NULL);
  }

  if (*nvar_chk_ec > 0) {
  dims_array[0] = *nvar_chk_ec * nx * (ny+k2d) * (nz+k3d);
  unk_e_x_tid = H5Tarray_create(H5T_NATIVE_FLOAT, rank_array, 
				 dims_array, NULL);
  dims_array[0] = *nvar_chk_ec * (nx+1) * ny * (nz+k3d);
  unk_e_y_tid = H5Tarray_create(H5T_NATIVE_FLOAT, rank_array, 
				 dims_array, NULL);
  dims_array[0] = *nvar_chk_ec * (nx+1) * (ny+k2d) * nz;
  unk_e_z_tid = H5Tarray_create(H5T_NATIVE_FLOAT, rank_array, 
				 dims_array, NULL);
  }

  if (*nvar_chk_nc > 0) {
  dims_array[0] = *nvar_chk_nc * (nx+1) * (ny+k2d) * (nz+k3d);
  unk_n_tid = H5Tarray_create(H5T_NATIVE_FLOAT, rank_array, 
				 dims_array, NULL);
  }

  offset = 0;

  H5Tinsert(*block_type_id, "lrefine",  
	    offset, H5T_NATIVE_FLOAT);
  offset = offset + sizeof(float);

  H5Tinsert(*block_type_id, "nodetype", 
	    offset, H5T_NATIVE_FLOAT);
  offset = offset + sizeof(float);

  H5Tinsert(*block_type_id, "which_child", 
	    offset, H5T_NATIVE_FLOAT);
  offset = offset + sizeof(float);

  H5Tinsert(*block_type_id, "gid", 
	    offset, gid_tid);
  offset = offset + (sizeof(float) * *ngid);

  H5Tinsert(*block_type_id, "bflags",  
	    offset, bflags_tid);
  offset = offset + (sizeof(float) * *mflags);

  H5Tinsert(*block_type_id, "coord", 
	    offset, coord_tid);
  offset = offset + (sizeof(float) * *mdim);

  H5Tinsert(*block_type_id, "bnd_box", 
	    offset , bnd_box_tid);
  offset = offset + (sizeof(float) * 2 * *mdim);

  H5Tinsert(*block_type_id, "work_block", 
	    offset, H5T_NATIVE_FLOAT);
  offset = offset + sizeof(float);

  if (*nvar_chk_cc > 0) {
    H5Tinsert(*block_type_id, "unk", 
	      offset, unk_tid);
    offset = offset + (sizeof(float) * *nvar_chk_cc * (*iu0-*il0) * (*ju0-*jl0) * (*ku0-*kl0));
  }
  if (*nvar_chk_fc > 0) {
    H5Tinsert(*block_type_id, "facevarx", 
	      offset, facevarx_tid);
    offset = offset + (sizeof(float) * *nvar_chk_fc * (*iu0-*il0+1) * (*ju0-*jl0) * (*ku0-*kl0));
    H5Tinsert(*block_type_id, "facevary", 
	      offset, facevary_tid);
    offset = offset + (sizeof(float) * *nvar_chk_fc * (*iu0-*il0) * (*ju0-*jl0+k2d) * (*ku0-*kl0));
    H5Tinsert(*block_type_id, "facevarz", 
	      offset, facevarz_tid);
    offset = offset + (sizeof(float) * *nvar_chk_fc * (*iu0-*il0) * (*ju0-*jl0) * (*ku0-*kl0+k3d));
  }
  if (*nvar_chk_ec > 0) {
    H5Tinsert(*block_type_id, "unk_e_x", 
	      offset, unk_e_x_tid);
    offset = offset + (sizeof(float) * *nvar_chk_ec * (*iu0-*il0) * (*ju0-*jl0+k2d) * (*ku0-*kl0+k3d));
    H5Tinsert(*block_type_id, "unk_e_y", 
	      offset, unk_e_y_tid);
    offset = offset + (sizeof(float) * *nvar_chk_ec * (*iu0-*il0+1) * (*ju0-*jl0) * (*ku0-*kl0+k3d));
    H5Tinsert(*block_type_id, "unk_e_z", 
	      offset, unk_e_z_tid);
    offset = offset + (sizeof(float) * *nvar_chk_ec * (*iu0-*il0+1) * (*ju0-*jl0+k2d) * (*ku0-*kl0));
  }
  if (*nvar_chk_nc > 0) {
    H5Tinsert(*block_type_id, "unk_n", 
	      offset, unk_n_tid);
  }

  H5Tclose(gid_tid);
  H5Tclose(bflags_tid);
  H5Tclose(coord_tid);
  H5Tclose(bnd_box_tid);
  if (*nvar_chk_cc > 0) {
    H5Tclose(unk_tid);
  }
  if (*nvar_chk_fc > 0) {
    H5Tclose(facevarx_tid);
    H5Tclose(facevary_tid);
    H5Tclose(facevarz_tid);
  }
  if (*nvar_chk_ec > 0) {
    H5Tclose(unk_e_x_tid);
    H5Tclose(unk_e_y_tid);
    H5Tclose(unk_e_z_tid);
  }
  if (*nvar_chk_nc > 0) {
    H5Tclose(unk_n_tid);
  }

}


void CREATE_Block_data_set_hdf5_r8(hid_t* dset_block, 
				  hid_t* file_id,
				  hid_t* block_type_id,
				  int*   tot_blocks_wr)
{

  hsize_t dims_block[1];
  hid_t   dataspace_block;

  /* Create the dataspace for block info */
  dims_block[0] = *tot_blocks_wr;
  dataspace_block = H5Screate_simple(1, dims_block, NULL);
  
  /* Create the dataset for block and cell centered data */
  /* Create a dataset collectively for each processor */
  *dset_block = H5Dcreate(*file_id, "Block Data", *block_type_id, 
                          dataspace_block, H5P_DEFAULT);
  H5Sclose(dataspace_block);

}

void block_pointer_to_var(float unk[], int iorder[], int nvar, int nx, int ny, int nz,
			  int kstart, int kend,
			  int jstart, int jend,
			  int istart, int iend,
			  int checkp_on[], int off1, int off2, int lb,
			  int *off, float *block_pointer)
{

  int udim[4]; int iii;
  for (iii = 0; iii < 4; iii++) {
    if (iorder[iii] == 1) {
      udim[iii] = nvar;
	  }
    if (iorder[iii] == 2) {
      udim[iii] = nx;
    }
    if (iorder[iii] == 3) {
      udim[iii] = ny;
    }
    if (iorder[iii] == 4) {
      udim[iii] = nz;
    }
  }
  
  int i; int j; int k; int ivar;
  for (k = kstart; k < kend; k++) {
    for (j = jstart; j < jend; j++) {
      for (i = istart; i < iend; i++) {
	for (ivar=0; ivar < nvar; ivar++) {
	  
	  int it[4];
	  for (iii = 0; iii < 4; iii++) {
	    if (iorder[iii] == 1) {
	      it[iii] = ivar;
	    }
	    if (iorder[iii] == 2) {
	      it[iii] = i;
	    }
	    if (iorder[iii] == 3) {
	      it[iii] = j;
	    }
	    if (iorder[iii] == 4) {
	      it[iii] = k;
	    }
	  }
	  
	  if (checkp_on[ivar + off1*ivar + off2] == 1) {
	    *(block_pointer+*off) =
	      unk[(udim[3] * udim[2] * udim[1] * udim[0] * lb) +
		  (udim[2] * udim[1] * udim[0] * it[3]) +