write_blocks_chombo_r4.c

做网格的好程序

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

#include "underscore.h"

void get_udim_chombo_r4(int nvar, int nx, int ny, int nz,
              int ivar, int istart, int jstart, int kstart,
              int iorder[4], int udim[4], int it[4]) {

  int iii;
  for (iii = 0; iii < 4; iii++) {
    if (iorder[iii] == 1) {
      udim[iii] = nvar;
      it[iii] = ivar;
    }
    if (iorder[iii] == 2) {
      udim[iii] = nx;
      it[iii] = istart;
    }
    if (iorder[iii] == 3) {
      udim[iii] = ny;
      it[iii] = jstart;
    }
    if (iorder[iii] == 4) {
      udim[iii] = nz;
      it[iii] = kstart;
    }
  }
}

void CREATE_File4 (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_Header4(hid_t*  file_id, 
		   int*    num_components,
		   int*    num_levels,
		   char*   compNames,
		   int*    ndim,
		   int*    nxb,
		   int*    nyb,
		   int*    nzb,
		   int     no_at_level[]) {

  hsize_t  dims[1];
  hid_t    attribute_id, dataspace_id, group_id, group_id2;
  hid_t    string_type, prob_domain_type_id;
  hid_t    dset_id;
  herr_t   status;
  int      mpi_rank;
  int      temp[1];
  float    rtemp[1];
  char     ctemp[20];
  char     name[20];
  int      i, ii;
  float    dx, dy, dz;

  typedef struct prob_domain_type_3d {
    int lo_i;
    int lo_j;
    int lo_k;
    int hi_i;
    int hi_j;
    int hi_k;
  } prob_domain_type_3d;
  prob_domain_type_3d prob_domain_3d[1];

  typedef struct prob_domain_type_2d {
    int lo_i;
    int lo_j;
    int hi_i;
    int hi_j;
  } prob_domain_type_2d;
  prob_domain_type_2d prob_domain_2d[1];

  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);

  group_id = H5Gopen(*file_id, "/");


  /* num_components */
  dataspace_id = H5Screate(H5S_SCALAR);
  attribute_id = H5Acreate(group_id, "num_components", 
			   H5T_NATIVE_INT, dataspace_id, 
			   H5P_DEFAULT);
  if (mpi_rank == 0) {
    status = H5Awrite(attribute_id, H5T_NATIVE_INT, num_components);
  }
  H5Sclose(dataspace_id);
  H5Aclose(attribute_id);


  /* iteration, this is a dummy */
  dataspace_id = H5Screate(H5S_SCALAR);
  attribute_id = H5Acreate(group_id, "iteration", 
			   H5T_NATIVE_INT, dataspace_id, 
			   H5P_DEFAULT);
  temp[0] = 0;
  if (mpi_rank == 0) {
    status = H5Awrite(attribute_id, H5T_NATIVE_INT, temp);
  }
  H5Sclose(dataspace_id);
  H5Aclose(attribute_id);


  /* num_levels */
  dataspace_id = H5Screate(H5S_SCALAR);
  attribute_id = H5Acreate(group_id, "num_levels", 
			   H5T_NATIVE_INT, dataspace_id, 
			   H5P_DEFAULT);
  if (mpi_rank == 0) {
    status = H5Awrite(attribute_id, H5T_NATIVE_INT, num_levels);
  }
  H5Sclose(dataspace_id);
  H5Aclose(attribute_id);


  /* time, this is a dummy */
  dataspace_id = H5Screate(H5S_SCALAR);
  attribute_id = H5Acreate(group_id, "time", 
			   H5T_NATIVE_FLOAT, dataspace_id, 
			   H5P_DEFAULT);
  rtemp[0] = 0;
  if (mpi_rank == 0) {
    status = H5Awrite(attribute_id, H5T_NATIVE_FLOAT, rtemp);
  }
  H5Sclose(dataspace_id);
  H5Aclose(attribute_id);


  /* compNames */
  for (i = 0; i < *num_components; i++) {
    dataspace_id = H5Screate(H5S_SCALAR);
    string_type  = H5Tcopy(H5T_C_S1);
    H5Tset_size(string_type, 20);
    sprintf(name, "component_%i",i);
    attribute_id = H5Acreate(group_id, name, 
			     string_type, dataspace_id, 
			     H5P_DEFAULT);
    for (ii = 0; ii < 20; ii++) {
      ctemp[ii] = compNames[i*20+ii];
    }
    if (mpi_rank == 0) {
      status = H5Awrite(attribute_id, string_type, ctemp);
    }
    H5Sclose(dataspace_id);
    H5Aclose(attribute_id);
    H5Tclose(string_type);
  }

  /* close group "/" */
  H5Gclose(group_id);

  /* Group Chombo global */
  group_id = H5Gcreate(*file_id, "/Chombo_global", 0);

  /* Add attribute "SpaceDim" to group "/Chombo_global" */
  dataspace_id = H5Screate(H5S_SCALAR);
  attribute_id = H5Acreate(group_id, "SpaceDim", 
			   H5T_NATIVE_INT, dataspace_id, 
			   H5P_DEFAULT);
  if (mpi_rank == 0) {
    status = H5Awrite(attribute_id, H5T_NATIVE_INT, ndim);
  }
  H5Sclose(dataspace_id);
  H5Aclose(attribute_id);

  /* Add attribute "testReal" to group "/Chombo_global" */
  /* This is a dummy */
  dataspace_id = H5Screate(H5S_SCALAR);
  attribute_id = H5Acreate(group_id, "testReal", 
			   H5T_NATIVE_FLOAT, dataspace_id, 
			   H5P_DEFAULT);
  rtemp[0] = 0;
  if (mpi_rank == 0) {
    status = H5Awrite(attribute_id, H5T_NATIVE_FLOAT, rtemp);
  }
  H5Sclose(dataspace_id);
  H5Aclose(attribute_id);

  /* close group "/Chombo_global" */
  H5Gclose(group_id);

  /* Create Groups for each level */
  for (i = 0; i < *num_levels; i++) {

    sprintf(name, "/level_%i",i);
    group_id = H5Gcreate(*file_id, name, 0);
    
    /* Add attribute for ref_ratio */
    dataspace_id = H5Screate(H5S_SCALAR);
    attribute_id = H5Acreate(group_id, "ref_ratio", 
			     H5T_NATIVE_INT, dataspace_id, 
			     H5P_DEFAULT);
    temp[0] = 2;
    if (mpi_rank == 0) {
      status = H5Awrite(attribute_id, H5T_NATIVE_INT, temp);
    }
    H5Sclose(dataspace_id);
    H5Aclose(attribute_id);
    
    /* Add attribute for "dt" */
    /* This is a dummy */
    dataspace_id = H5Screate(H5S_SCALAR);
    attribute_id = H5Acreate(group_id, "dt", 
			     H5T_NATIVE_FLOAT, dataspace_id, 
			     H5P_DEFAULT);
    rtemp[0] = 0.;
    if (mpi_rank == 0) {
      status = H5Awrite(attribute_id, H5T_NATIVE_FLOAT, rtemp);
    }
    H5Sclose(dataspace_id);
    H5Aclose(attribute_id);
    
    /* Add attribute for "dx" */
    dataspace_id = H5Screate(H5S_SCALAR);
    attribute_id = H5Acreate(group_id, "dx", 
			     H5T_NATIVE_FLOAT, dataspace_id, 
			     H5P_DEFAULT);
    /* This is not really right, but things should look OK 
       as long as the cell sizes change by only a factor of 2
       from one level to the next
       This will change the total domain size 
       What we really want is a 1 x 1 x 1 box in all cases, but relative
       sizes should be OK .
       Maybe its best just to assume that the top level is just 1 block ??? */
    dx = 1./(pow(2,i)*(float) *nxb);
    if (mpi_rank == 0) {
      status = H5Awrite(attribute_id, H5T_NATIVE_FLOAT, &dx);
    }
    H5Sclose(dataspace_id);
    H5Aclose(attribute_id);

    /* Add attribute for "prob_domain" */

    if (*ndim == 3) {

    prob_domain_type_id = H5Tcreate (H5T_COMPOUND, sizeof(prob_domain_type_3d));

    H5Tinsert(prob_domain_type_id, "lo_i",  
	      HOFFSET(prob_domain_type_3d, lo_i), H5T_NATIVE_INT);
    H5Tinsert(prob_domain_type_id, "lo_j",  
	      HOFFSET(prob_domain_type_3d, lo_j), H5T_NATIVE_INT);
    H5Tinsert(prob_domain_type_id, "lo_k",  
	      HOFFSET(prob_domain_type_3d, lo_k), H5T_NATIVE_INT);
    H5Tinsert(prob_domain_type_id, "hi_i",  
	      HOFFSET(prob_domain_type_3d, hi_i), H5T_NATIVE_INT);
    H5Tinsert(prob_domain_type_id, "hi_j",  
	      HOFFSET(prob_domain_type_3d, hi_j), H5T_NATIVE_INT);
    H5Tinsert(prob_domain_type_id, "hi_k",  
	      HOFFSET(prob_domain_type_3d, hi_k), H5T_NATIVE_INT);

    dims[0] = 1;
    dataspace_id = H5Screate(H5S_SCALAR);
    attribute_id = H5Acreate(group_id, "prob_domain", 
			     prob_domain_type_id, dataspace_id, 
			     H5P_DEFAULT);

    prob_domain_3d[0].lo_i = 0;
    prob_domain_3d[0].lo_j = 0;
    prob_domain_3d[0].lo_k = 0;
    prob_domain_3d[0].hi_i = (pow(2,i) * *nxb) - 1;
    prob_domain_3d[0].hi_j = (pow(2,i) * *nyb) - 1;
    prob_domain_3d[0].hi_k = (pow(2,i) * *nzb) - 1;
    
    if (mpi_rank == 0) {
      status = H5Awrite(attribute_id, prob_domain_type_id, prob_domain_3d);
    }
    H5Sclose(dataspace_id);
    H5Aclose(attribute_id);

    } else {

    prob_domain_type_id = H5Tcreate (H5T_COMPOUND, sizeof(prob_domain_type_2d));

    H5Tinsert(prob_domain_type_id, "lo_i",  
	      HOFFSET(prob_domain_type_2d, lo_i), H5T_NATIVE_INT);
    H5Tinsert(prob_domain_type_id, "lo_j",  
	      HOFFSET(prob_domain_type_2d, lo_j), H5T_NATIVE_INT);
    H5Tinsert(prob_domain_type_id, "hi_i",  
	      HOFFSET(prob_domain_type_2d, hi_i), H5T_NATIVE_INT);
    H5Tinsert(prob_domain_type_id, "hi_j",  
	      HOFFSET(prob_domain_type_2d, hi_j), H5T_NATIVE_INT);

    dims[0] = 1;
    dataspace_id = H5Screate(H5S_SCALAR);
    attribute_id = H5Acreate(group_id, "prob_domain", 
			     prob_domain_type_id, dataspace_id, 
			     H5P_DEFAULT);

    prob_domain_2d[0].lo_i = 0;
    prob_domain_2d[0].lo_j = 0;
    prob_domain_2d[0].hi_i = (pow(2,i) * *nxb) - 1;
    prob_domain_2d[0].hi_j = (pow(2,i) * *nyb) - 1;
    
    if (mpi_rank == 0) {
      status = H5Awrite(attribute_id, prob_domain_type_id, prob_domain_2d);
    }
    H5Sclose(dataspace_id);
    H5Aclose(attribute_id);

    }

    /* Add "boxes" dataset */
    dims[0] = no_at_level[i];
    dataspace_id = H5Screate_simple(1, dims, NULL);
    dset_id = H5Dcreate(group_id, "boxes", prob_domain_type_id, 
			dataspace_id, H5P_DEFAULT);
    H5Dclose(dset_id);
    H5Sclose(dataspace_id);

    H5Tclose(prob_domain_type_id);

    /* Add "data:datatype=0" dataset */
    dims[0] = *num_components * *nxb * *nyb * *nzb * no_at_level[i];
    dataspace_id = H5Screate_simple(1, dims, NULL);
    dset_id = H5Dcreate(group_id, "data:datatype=0", H5T_NATIVE_FLOAT, 
			dataspace_id, H5P_DEFAULT);
    H5Dclose(dset_id);
    H5Sclose(dataspace_id);

    /* Group "data_attributes" */
    group_id2 = H5Gcreate(group_id, "data_attributes", 0);
    
    /* Add attribute "comps" to group "data_attributes" */
    dataspace_id = H5Screate(H5S_SCALAR);
    attribute_id = H5Acreate(group_id2, "comps", 
			     H5T_NATIVE_INT, dataspace_id, 
			     H5P_DEFAULT);
    if (mpi_rank == 0) {
      status = H5Awrite(attribute_id, H5T_NATIVE_INT, num_components);
    }
    H5Sclose(dataspace_id);
    H5Aclose(attribute_id);

    /* close group "data_attributes" */
    H5Gclose(group_id2);

    /* close the group for level i */
    H5Gclose(group_id);

  } /* end loop over levels */

}


void WRITE_Blocks4(hid_t* file_id, 
		   int* num_levels,
		   int* num_components,
		   int* nxb,
		   int* nyb,
		   int* nzb,
		   int* max_lnblocks, 
		   int* n_to_left, 
		   int  n_to_left_level[], 
		   int* lnblocks, 
		   int* mdim,
		   int* ndim,
		   int* mflags,
		   int lrefine[], 
		   int nodetype[], 
		   int which_child[], 
		   int bflags[],
		   float coord[],
		   float bnd_box[],
		   float work_block[],
		   float unk[],
		   int*  nvar,
		   int*  nvar_chk_cc,
		   int   checkp_on_cc[],
		   float facevarx[], float facevary[], float facevarz[],
		   int*  nbndvar,
		   int*  nvar_chk_fc,
		   int   checkp_on_fc[],
		   float unk_e_x[], float unk_e_y[], float unk_e_z[],
		   int*  nbndvare,
		   int*  nvar_chk_ec,
		   int   checkp_on_ec[],
		   float unk_n[],
		   int*  nbndvarc,
		   int*  nvar_chk_nc,
		   int   checkp_on_nc[],
		   int*  nx,  int* ny, int* nz,
		   int*  il0, int* iu0,
		   int*  jl0, int* ju0, 
		   int*  kl0, int* ku0,
		   int   iorder[])
{
  hid_t  memspace, filespace;
  hsize_t  count[1];
  int rank_array;
  hsize_t  dims[1], dims_array[1];
  hssize_t offset[1];
  hid_t plist_id;
  hid_t dset_id, dataspace_id, group_id, prob_domain_type_id;
  int ivar, i, j, k, lb, i2, j2, k2, nx2, ny2, ivar2;
  int istart, iend, jstart, jend, kstart, kend, k3d;
  herr_t status;
  char name[20];
  int  n_at_level[*num_levels];
  int  block_no_at_level[*lnblocks];
  float xmin, xmax, ymin, ymax, zmin, zmax, dx;
  int lb2;
  float unk_t[*nxb * *nyb * *nzb];
  int mpi_rank;
  int ii, jj, kk;
  float xl, yl, zl;
  float xleft, yleft, zleft;


  typedef struct prob_domain_type_3d {