1 /*   Program read_grid_str_parinzone.c    */
2 /*
3 Reads a simple 3-D structured grid from a CGNS file.  Each processor reads a
4 slab of data from one zone (parallelism within a zone).
5 
6 The CGNS grid file 'grid_piz_c.cgns' must already exist.
7 
8 mpicxx read_grid_str_parinzone.c -lcgns -lhdf5 -lsz -lz -o read_grid_str_parinzone
9 mpirun -np 2 write_grid_str_parinzone
10 */
11 
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <string.h>
15 
16 #include "pcgnslib.h"
17 #include "mpi.h"
18 
main(int argc,const char * argv[])19 int main(int argc, const char* argv[])
20 {
21    cgsize_t zoneSize[3][3];
22    int n, comm_size, comm_rank;
23    int index_file, index_base, index_zone;
24    char zonename[33];
25 
26    MPI_Init(&argc, (char***)(&argv));
27    MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
28    MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
29    cgp_mpi_comm(MPI_COMM_WORLD);
30 
31 /* open CGNS file for read-only */
32    if (cgp_open("grid_piz_c.cgns", CG_MODE_READ, &index_file)) cg_error_exit();
33 /* we know there is only one base (real working code would check!) */
34     index_base = 1;
35 /* we know there is only one zone (real working code would check!) */
36     index_zone = 1;
37 /* get zone size (and name - although not needed here) */
38     if (cg_zone_read(index_file, index_base, index_zone, zonename,
39                      (cgsize_t*)zoneSize)) cg_error_exit();
40 
41 /* COLLECTIVE READING OF FILE DATA -- each processor reads a part of the zone */
42 
43    /* partition the k-index of the zone among the processes */
44    int kIdxBeg;  /* beginning k-index for this processor */
45    int numLocalkIdx = zoneSize[0][2]/comm_size;  /* the number of k-indices
46                                                     local to this process */
47    {
48      int numUnevenkIdx = zoneSize[0][2] - numLocalkIdx*comm_size;
49      if (comm_rank < numUnevenkIdx)
50        {
51          ++numLocalkIdx;
52          kIdxBeg = comm_rank*numLocalkIdx;
53        }
54      else
55        {
56          kIdxBeg = numUnevenkIdx*(numLocalkIdx + 1);
57          kIdxBeg += (comm_rank - numUnevenkIdx)*numLocalkIdx;
58        }
59    }
60 
61    int s_rmin[3], s_rmax[3], m_dimvals[3], m_rmin[3], m_rmax[3];
62    double *x = NULL;
63    double *y = NULL;
64    double *z = NULL;
65    int ni, nj, nk;
66    if (numLocalkIdx > 0)
67      {
68        /* allocate memory for reading */
69        ni = zoneSize[0][0];
70        nj = zoneSize[0][1];
71        nk = numLocalkIdx;
72        const int num_vertex = ni*nj*nk;
73        x = (double*)malloc(3*num_vertex*sizeof(double));
74        y = x + num_vertex;
75        z = y + num_vertex;
76        // shape in file space
77        for (n = 0; n < 2; ++n)
78          {
79            s_rmin[n] = 1;
80            s_rmax[n] = zoneSize[0][n];
81          }
82        s_rmin[2] = kIdxBeg + 1;
83        s_rmax[2] = kIdxBeg + numLocalkIdx;
84        // shape in memory
85        for (n = 0; n < 2; ++n)
86          {
87            m_dimvals[n] = zoneSize[0][n];
88            m_rmin[n]    = 1;
89            m_rmax[n]    = zoneSize[0][n];
90          }
91        m_dimvals[2] = numLocalkIdx;
92        m_rmin[2]    = 1;
93        m_rmax[2]    = numLocalkIdx;
94      }
95    /* the file data is defined as single but will be read as double */
96    if (cgp_coord_general_read_data(index_file, index_base, index_zone, 1,
97                                    s_rmin, s_rmax, CGNS_ENUMV(RealDouble),
98                                    3, m_dimvals, m_rmin, m_rmax,
99                                    x)) cgp_error_exit();
100    if (cgp_coord_general_read_data(index_file, index_base, index_zone, 2,
101                                    s_rmin, s_rmax, CGNS_ENUMV(RealDouble),
102                                    3, m_dimvals, m_rmin, m_rmax,
103                                    y)) cgp_error_exit();
104    if (cgp_coord_general_read_data(index_file, index_base, index_zone, 3,
105                                    s_rmin, s_rmax, CGNS_ENUMV(RealDouble),
106                                    3, m_dimvals, m_rmin, m_rmax,
107                                    z)) cgp_error_exit();
108    if (8 >= kIdxBeg && 8 < kIdxBeg + numLocalkIdx)
109      {
110        printf("\nSuccessfully read grid from file grid_piz_c.cgns\n");
111        printf("  For example, from process %d, zone 1 x,y,z[8][18][16]= "
112               "%f, %f, %f\n", comm_rank,
113               x[((8 - kIdxBeg)*nj + 18)*ni + 16],
114               y[((8 - kIdxBeg)*nj + 18)*ni + 16],
115               z[((8 - kIdxBeg)*nj + 18)*ni + 16]);
116      }
117    if (numLocalkIdx > 0)
118      {
119        free(x);
120      }
121 /* close CGNS file */
122    cgp_close(index_file);
123    MPI_Finalize();
124    if (comm_rank == 0)
125      {
126        printf("\nProgram successful... ending now\n");
127      }
128    return 0;
129 }
130