c: simple read of gbpCode grid

read_grid.c

#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <math.h>

#define SNAPSHOT 99
#define INPUT_DIR "/home/gpoole/SSimPL/simulations/GiggleZ/GiggleZ_HR/grids/"

/**
 * Compile with:
 *  gcc -std=c99 -O2 -Wall read_grid.c -o read_grid
 **/ 

static inline int index(int i, int j, int k, int dim, int type)
{
  int ind;
  // type =0 -> real
  // type =1 -> padded
  switch(type)
  {
    case 1:
      ind  = k + (2*(dim/2 +1)) * (j + dim * (i));
      break;
    case 0:
      ind = k + dim * (j + dim *(i));
      break;
  }

  return ind;
}


static inline void read_identifier(FILE *fin, bool skip_flag)
{
    char identifier[32];

    fread(identifier, sizeof(identifier), 1, fin);

    if (skip_flag)
        fprintf(stderr, "Skipping grid: %s...\n", identifier);
    else
        printf("Reading grid: %s...\n", identifier);
}


static int read_grid(
    int            i_grid,
    int           *dim,
    float        **grid)
{
    char   fname[512];
    FILE   *fin;
    int    n_cell[3];
    double box_size[3];
    int    n_grids;
    int    ma_scheme;
    int    n_elem;

    // Construct the input filename
    sprintf(fname, "%s/snapshot_%03d_dark_grid.dat", INPUT_DIR, SNAPSHOT);
    // ... and open
    fin = fopen(fname, "rb");
    if (!fin)
    {
        fprintf(stderr, "Failed to open file: %s", fname);
        return(EXIT_FAILURE);
    }

    // Read the header
    fread(n_cell, sizeof(int), 3, fin);
    fread(box_size, sizeof(double), 3, fin);
    fread(&n_grids, sizeof(int), 1, fin);
    fread(&ma_scheme, sizeof(int), 1, fin);

    printf("Reading grid for snapshot %d\n", SNAPSHOT);
    printf("n_cell = [%d, %d, %d]\n", n_cell[0], n_cell[1], n_cell[2]);
    printf("box_size = [%.2f, %.2f, %.2f]\n", box_size[0], box_size[1], box_size[2]);
    printf("ma_scheme = %d\n", ma_scheme);


    if (n_grids != 4)
    {
        fprintf(stderr, "n_grids != 4 as expected...\n");
        fclose(fin);
        return -1;
    }

    *dim = n_cell[0];
    n_elem = n_cell[0] * n_cell[1] * n_cell[2];

    // Read the grid
    for (int ii = 0; ii < i_grid; ii++)
    {
        read_identifier(fin, true);
        fseek(fin, sizeof(float) * n_elem, SEEK_CUR);
    }
    read_identifier(fin, false);

    // Malloc the grid and read
    *grid = calloc(n_elem, sizeof(double));
    for (int ii = 0; ii < n_elem; ii++)
      fread(&((*grid)[ii]), sizeof(float), 1, fin);

    // Close the file
    fclose(fin);

    return 0;
}


int main(int argc, char *argv[])
{
  int dim;
  float *grid;
  int test_ind[3] = {0, 255, 511};

  // read the grid (0 is the density grid)
  read_grid(0, &dim, &grid);

  // print a few values
  printf("\n");
  printf("-----------\n");
  
  for (int ii = 0; ii < 3; ii++)
  {
    int jj = test_ind[ii];
    printf("grid[%03d,%03d,%03d]\t= %g\n", jj, jj, jj, grid[index(jj,jj,jj,dim,0)]);
  }
  
  printf("-----------\n");

  free(grid);

  return 0;
}