write_blocks_chombo_r8.c

做网格的好程序

  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 boxes_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 boxes_2d[1];
  

  /* Create property list for collective dataset write */
  plist_id = H5Pcreate(H5P_DATASET_XFER);
  /*H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);*/
  H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT);

  /* Create "prob_domain_id" */

  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);

  } 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);

  }
  
  nx2 = *iu0 - *il0;
  ny2 = *ju0 - *jl0;
  kstart = *kl0;
  kend   = *ku0;
  jstart = *jl0;
  jend   = *ju0;
  istart = *il0;
  iend   = *iu0;
  k3d    = 1;
  if (*ndim < 3) {
    kstart = 0;
    kend   = 1;
    k3d    = 0;
  }
  if (*ndim < 2) {
    kstart = 0;
    kend   = 1;
    jstart = 0;
    jend   = 1;
    k3d    = 0;
  }
  
  /* compute no of blocks at each level */
  for (i = 0; i < *num_levels; i++) {
    n_at_level[i] = 0;
  }
  for (i = 0; i < *lnblocks; i++) {
    block_no_at_level[i] = 0;
  }

  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
  for (lb = 0; lb < *lnblocks; lb++) {
    block_no_at_level[lb] = n_to_left_level[lrefine[lb]-1] + 
                            n_at_level[lrefine[lb]-1];
    n_at_level[lrefine[lb]-1] = n_at_level[lrefine[lb]-1] + 1;
  }

  /* Find the minimum (left) sides of the domain */

  xl = 1.e30;
  yl = 1.e30;
  zl = 1.e30;
  for (lb = 0; lb < *lnblocks; lb++) {
    if (bnd_box[lb* *mdim*2 + 2*0 + 0] < xl) {
      xl = bnd_box[lb* *mdim*2 + 2*0 + 0];
    }
    if (bnd_box[lb* *mdim*2 + 2*1 + 0] < yl) {
      yl = bnd_box[lb* *mdim*2 + 2*1 + 0];
    }
    if (bnd_box[lb* *mdim*2 + 2*2 + 0] < zl) {
      zl = bnd_box[lb* *mdim*2 + 2*2 + 0];
    }
  }

  status = MPI_Allreduce(&xl, &xleft, 1, MPI_FLOAT, MPI_MIN, MPI_COMM_WORLD);
  status = MPI_Allreduce(&yl, &yleft, 1, MPI_FLOAT, MPI_MIN, MPI_COMM_WORLD);
  status = MPI_Allreduce(&zl, &zleft, 1, MPI_FLOAT, MPI_MIN, MPI_COMM_WORLD);

  for (lb = 0; lb < *max_lnblocks; lb++) { 

    if (lb < *lnblocks) {
      lb2 = lb;
    } else {
      lb2 = *lnblocks - 1;
    }

    /* Open Group for this level */
    sprintf(name, "/level_%i",lrefine[lb2]-1);
    group_id = H5Gopen(*file_id, name);


    /* Open "boxes" dataset */
    dims[0] = 1;
    dataspace_id = H5Screate_simple(1, dims, NULL);
    dset_id = H5Dopen(group_id, "boxes");
    
    offset[0] = block_no_at_level[lb2]; 
    count[0] = 1;
    filespace = H5Dget_space(dset_id);
    H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, 
			NULL, count, NULL);
    
    dx      = bnd_box[lb2* *mdim*2 + 2*0 + 1] - bnd_box[lb2* *mdim*2 + 2*0 + 0];
    dx      = dx/(float) *nxb;
    xmin    = bnd_box[lb2* *mdim*2 + 2*0 + 0] + dx/2.;
    xmax    = bnd_box[lb2* *mdim*2 + 2*0 + 1] - dx/2.;
    ymin    = bnd_box[lb2* *mdim*2 + 2*1 + 0] + dx/2.;
    ymax    = bnd_box[lb2* *mdim*2 + 2*1 + 1] - dx/2.;
    zmin    = bnd_box[lb2* *mdim*2 + 2*2 + 0] + dx/2.;
    zmax    = bnd_box[lb2* *mdim*2 + 2*2 + 1] - dx/2.;

    if (*ndim == 3) {

    boxes_3d[0].lo_i = ((int)((xmin - xleft)/dx));
    boxes_3d[0].lo_j = ((int)((ymin - yleft)/dx));
    boxes_3d[0].lo_k = ((int)((zmin - zleft)/dx));
    boxes_3d[0].hi_i = ((int)((xmax - xleft)/dx));
    boxes_3d[0].hi_j = ((int)((ymax - yleft)/dx));
    boxes_3d[0].hi_k = ((int)((zmax - zleft)/dx));

    status = H5Dwrite(dset_id, prob_domain_type_id, dataspace_id, filespace,
		      plist_id, boxes_3d);

    } else {

    boxes_2d[0].lo_i = ((int)((xmin - xleft)/dx));
    boxes_2d[0].lo_j = ((int)((ymin - yleft)/dx));
    boxes_2d[0].hi_i = ((int)((xmax - xleft)/dx));
    boxes_2d[0].hi_j = ((int)((ymax - yleft)/dx));

    status = H5Dwrite(dset_id, prob_domain_type_id, dataspace_id, filespace,
		      plist_id, boxes_2d);

    }

    /* close group data space and data */
    H5Sclose(dataspace_id);
    H5Dclose(dset_id);
    H5Sclose(filespace);

    status = MPI_Barrier(MPI_COMM_WORLD);

    /*****************************************************************/

    /* Cell Centered Data */
    
    count[0]   = *nxb * *nyb * *nzb;
    offset[0]  = block_no_at_level[lb2] * count[0] * *num_components; 

    if (*nvar_chk_cc > 0) {
      
      ivar2 = 0;
      for (ivar=0; ivar < *nvar; ivar++) {

	if (checkp_on_cc[ivar] == 1) {

	  /* Open "data:datatype=0" dataset */
	  dims[0] = *nxb * *nyb * *nzb;
	  dataspace_id = H5Screate_simple(1, dims, NULL);
	  dset_id = H5Dopen(group_id, "data:datatype=0");
	  
	  filespace = H5Dget_space(dset_id);
	  H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, 
			      NULL, count, NULL);

	  offset[0] = offset[0] + count[0];
    
	  k2 = 0;
	  for (k = kstart; k < kend; k++) {
	    j2 = 0;
	    for (j = jstart; j < jend; j++) {
	      i2 = 0;
	      for (i = istart; i < iend; i++) {
		
		int udim[4]; int it[4];
		get_udim_chombo (*nvar, *nx, *ny, *nz, ivar, i, j, k, iorder, udim, it);

		unk_t[i2 + 
		      nx2 * j2 + 
		      nx2 * ny2 * k2] =
		  unk[(udim[3] * udim[2] * udim[1] * udim[0] * lb2) + 
		      (udim[2] * udim[1] * udim[0] * it[3]) + 
		      (udim[1] * udim[0] * it[2]) + 
		      (udim[0] * it[1]) + 
		      it[0]];

		i2 = i2 + 1;
	      } /* end i loop */
	      j2 = j2 + 1;
	    } /* end j loop */
	    k2 = k2 + 1;
	  } /* end k loop */


	  status = H5Dwrite(dset_id, H5T_NATIVE_FLOAT, dataspace_id, filespace, 
			    plist_id, unk_t);

	  H5Sclose(dataspace_id);
	  H5Dclose(dset_id);
	  H5Sclose(filespace);

	  ivar2 = ivar2 + 1;
	} /* end if (checkp */
      } /* end ivar loop */

    } /* end if (nvar_chk_cc */

    /* Face Centered Data */


    if (*nvar_chk_fc > 0) {
      
      /* facevarx */

      ivar2 = 0;
      for (ivar=0; ivar < *nbndvar; ivar++) {
	if (checkp_on_fc[ivar] == 1) {

	  /* Open "data:datatype=0" dataset */
	  dims[0] = *nxb * *nyb * *nzb;
	  dataspace_id = H5Screate_simple(1, dims, NULL);
	  dset_id = H5Dopen(group_id, "data:datatype=0");
	  
	  filespace = H5Dget_space(dset_id);
	  H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, 
			      NULL, count, NULL);

	  offset[0] = offset[0] + count[0]; 
    
	  k2 = 0;
	  for (k = kstart; k < kend; k++) {
	    j2 = 0;
	    for (j = jstart; j < jend; j++) {
	      i2 = 0;
	      for (i = istart; i < iend; i++) {
		
		int udim[4]; int it[4]; int it2[4];
		get_udim_chombo (*nbndvar, *nx+1, *ny, *nz, ivar, i, j, k, iorder, udim, it);
		get_udim_chombo (*nbndvar, *nx+1, *ny, *nz, ivar, i+1, j, k, iorder, udim, it2);

		unk_t[i2 + 
		      nx2 * j2 + 
		      nx2 * ny2 * k2] =
		  facevarx[(udim[3] * udim[2] * udim[1] * udim[0] * lb2) + 
			   (udim[2] * udim[1] * udim[0] * it[3]) + 
			   (udim[1] * udim[0] * it[2]) + 
			   (udim[0] * it[1]) + 
			   it[0]] +
		  facevarx[(udim[3] * udim[2] * udim[1] * udim[0] * lb2) + 
			   (udim[2] * udim[1] * udim[0] * it2[3]) + 
			   (udim[1] * udim[0] * it2[2]) + 
			   (udim[0] * it2[1]) + 
			   it2[0]];

		unk_t[i2 + 
		      nx2 * j2 + 
		      nx2 * ny2 * k2] =
		  unk_t[i2 + 
			nx2 * j2 + 
			nx2 * ny2 * k2] / 2.;

		i2 = i2 + 1;
	      } /* end i loop */
	      j2 = j2 + 1;
	    } /* end j loop */
	    k2 = k2 + 1;
	  } /* end k loop */

	  status = H5Dwrite(dset_id, H5T_NATIVE_FLOAT, dataspace_id, filespace, 
			    plist_id, unk_t);

	  H5Sclose(dataspace_id);
	  H5Dclose(dset_id);
	  H5Sclose(filespace);

	  ivar2 = ivar2 + 1;
	} /* end if (checkp */
      } /* end ivar loop */

      /* facevary */

      ivar2 = 0;
      for (ivar=0; ivar < *nbndvar; ivar++) {
	if (checkp_on_fc[ivar] == 1) {

	  /* Open "data:datatype=0" dataset */
	  dims[0] = *nxb * *nyb * *nzb;
	  dataspace_id = H5Screate_simple(1, dims, NULL);
	  dset_id = H5Dopen(group_id, "data:datatype=0");
	  
	  filespace = H5Dget_space(dset_id);
	  H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, 
			      NULL, count, NULL);

	  offset[0] = offset[0] + count[0]; 
    
	  k2 = 0;
	  for (k = kstart; k < kend; k++) {
	    j2 = 0;
	    for (j = jstart; j < jend; j++) {
	      i2 = 0;
	      for (i = istart; i < iend; i++) {
		
		int udim[4]; int it[4]; int it2[4];
		get_udim_chombo (*nbndvar, *nx, *ny+1, *nz, ivar, i, j, k, iorder, udim, it);
		get_udim_chombo (*nbndvar, *nx, *ny+1, *nz, ivar, i+1, j+1, k, iorder, udim, it2);

		unk_t[i2 + 
		      nx2 * j2 + 
		      nx2 * ny2 * k2] =
		  facevary[(udim[3] * udim[2] * udim[1] * udim[0] * lb2) + 
			   (udim[2] * udim[1] * udim[0] * it[3]) + 
			   (udim[1] * udim[0] * it[2]) + 
			   (udim[0] * it[1]) + 
			   it[0]] +
		  facevary[(udim[3] * udim[2] * udim[1] * udim[0] * lb2) + 
			   (udim[2] * udim[1] * udim[0] * it2[3]) + 
			   (udim[1] * udim[0] * it2[2]) + 
			   (udim[0] * it2[1]) + 
			   it2[0]];

		unk_t[i2 + 
		      nx2 * j2 + 
		      nx2 * ny2 * k2] =
		  unk_t[i2 + 
			nx2 * j2 + 
			nx2 * ny2 * k2] / 2.;

		i2 = i2 + 1;
	      } /* end i loop */
	      j2 = j2 + 1;
	    } /* end j loop */
	    k2 = k2 + 1;
	  } /* end k loop */

	  status = H5Dwrite(dset_id, H5T_NATIVE_FLOAT, dataspace_id, filespace, 
			    plist_id, unk_t);

	  H5Sclose(dataspace_id);
	  H5Dclose(dset_id);
	  H5Sclose(filespace);

	  ivar2 = ivar2 + 1;
	} /* end if (checkp */
      } /* end ivar loop */

      /* facevarz */

      ivar2 = 0;
      for (ivar=0; ivar < *nbndvar; ivar++) {
	if (checkp_on_fc[ivar] == 1) {

	  /* Open "data:datatype=0" dataset */
	  dims[0] = *nxb * *nyb * *nzb;
	  dataspace_id = H5Screate_simple(1, dims, NULL);
	  dset_id = H5Dopen(group_id, "data:datatype=0");
	  
	  filespace = H5Dget_space(dset_id);
	  H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, 
			      NULL, count, NULL);

	  offset[0] = offset[0] + count[0]; 
    
	  k2 = 0;
	  for (k = kstart; k < kend; k++) {
	    j2 = 0;
	    for (j = jstart; j < jend; j++) {
	      i2 = 0;
	      for (i = istart; i < iend; i++) {
		
		int udim[4]; int it[4]; int it2[4];
		get_udim_chombo (*nbndvar, *nx, *ny, *nz+k3d, ivar, i, j, k, iorder, udim, it);
		get_udim_chombo (*nbndvar, *nx, *ny, *nz+k3d, ivar, i, j, k+k3d, iorder, udim, it2);

		unk_t[i2 + 
		      nx2 * j2 + 
		      nx2 * ny2 * k2] =
		  facevarz[(udim[3] * udim[2] * udim[1] * udim[0] * lb2) + 
			   (udim[2] * udim[1] * udim[0] * it[3]) + 
			   (udim[1] * udim[0] * it[2]) + 
			   (udim[0] * it[1]) + 
			   it[0]] +
		  facevarz[(udim[3] * udim[2] * udim[1] * udim[0] * lb2) + 
			   (udim[2] * udim[1] * udim[0] * it2[3]) + 
			   (udim[1] * udim[0] * it2[2]) + 
			   (udim[0] * it2[1]) + 
			   it2[0]];

		unk_t[i2 + 
		      nx2 * j2 + 
		      nx2 * ny2 * k2] =
		  unk_t[i2 + 
			nx2 * j2 + 
			nx2 * ny2 * k2] / 2.;

		i2 = i2 + 1;
	      } /* end i loop */
	      j2 = j2 + 1;
	    } /* end j loop */
	    k2 = k2 + 1;
	  } /* end k loop */

	  status = H5Dwrite(dset_id, H5T_NATIVE_FLOAT, dataspace_id, filespace, 
			    plist_id, unk_t);

	  H5Sclose(dataspace_id);
	  H5Dclose(dset_id);
	  H5Sclose(filespace);