ReadArraywInterp.c

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <mpivars.h>
#include <hypar.h>

static int ReadArraywInterpSerial(int,int,int*,int*,int*,int,void*,void*,double*,double*,char*,int*);

int ReadArraywInterp( int     ndims,          
                      int     nvars,          
                      int     *dim_global,    
                      int     *dim_local,     
                      int     *dim_global_src,
                      int     ghosts,         
                      void    *s,             
                      void    *m,             
                      double  *x,             
                      double  *u,             
                      char    *fname_root,    
                      int     *read_flag      
                     )
{
  HyPar         *solver = (HyPar*) s;
  MPIVariables  *mpi    = (MPIVariables*) m;
  _DECLARE_IERR_;

  int retval = ReadArraywInterpSerial(  ndims,
                                        nvars,
                                        dim_global,
                                        dim_local,
                                        dim_global_src,
                                        ghosts,
                                        s,
                                        m,
                                        x,
                                        u,
                                        fname_root,
                                        read_flag );
  if (retval) return retval;

  if (x) {
    int offset, d;
    /* exchange MPI-boundary values of x between processors */
    offset = 0;
    for (d = 0; d < ndims; d++) {
      IERR MPIExchangeBoundaries1D(mpi,&x[offset],dim_local[d],
                                   ghosts,d,ndims); CHECKERR(ierr);
      offset  += (dim_local [d] + 2*ghosts);
    }
    /* fill in ghost values of x at physical boundaries by extrapolation */
    offset = 0;
    for (d = 0; d < ndims; d++) {
      double *X     = &x[offset];
      int    *dim   = dim_local, i;
      if (mpi->ip[d] == 0) {
        /* fill left boundary along this dimension */
        for (i = 0; i < ghosts; i++) {
          int delta = ghosts - i;
          X[i] = X[ghosts] + ((double) delta) * (X[ghosts]-X[ghosts+1]);
        }
      }
      if (mpi->ip[d] == mpi->iproc[d]-1) {
        /* fill right boundary along this dimension */
        for (i = dim[d]+ghosts; i < dim[d]+2*ghosts; i++) {
          int delta = i - (dim[d]+ghosts-1);
          X[i] =  X[dim[d]+ghosts-1] 
                  + ((double) delta) * (X[dim[d]+ghosts-1]-X[dim[d]+ghosts-2]);
        }
      }
      offset  += (dim[d] + 2*ghosts);
    }
  }

  return 0;
}

int ReadArraywInterpSerial( int     ndims,          
                            int     nvars,          
                            int     *dim_global,    
                            int     *dim_local,     
                            int     *dim_global_src,
                            int     ghosts,         
                            void    *s,             
                            void    *m,             
                            double  *x,             
                            double  *u,             
                            char    *fname_root,    
                            int     *read_flag      
                          )
{
  HyPar         *solver = (HyPar*)        s;
  MPIVariables  *mpi    = (MPIVariables*) m;
  int           i, d, ferr, index[ndims];
  double        *ug_src = NULL, *xg_src = NULL;
  double        *ug = NULL, *xg = NULL;
  _DECLARE_IERR_;

  *read_flag = 0;
  /* Only root process reads from the file */
  if (!mpi->rank) {

    /* read data from file - this data is of dimensions given by dim_global_sec */
    if (!strcmp(solver->ip_file_type,"ascii")) {
      char filename[_MAX_STRING_SIZE_];
      strcpy(filename,fname_root);
      strcat(filename,".inp");
      FILE *in; in = fopen(filename,"r");
      if (!in) *read_flag = 0;
      else {
        *read_flag = 1;
        /* Reading from file */
        printf("Reading array from ASCII file %s (Serial mode).\n",filename);
        int size,offset;
        /* allocate global solution array */
        size   = 1; for (d=0; d<ndims; d++) size *= dim_global_src[d]; size *= nvars;
        ug_src = (double*) calloc(size,sizeof(double));
        size   = 0; for (d=0; d<ndims; d++) size += dim_global_src[d];
        xg_src = (double*) calloc(size,sizeof(double));

        /* read grid */
        offset = 0;
        for (d = 0; d < ndims; d++) {
          for (i = 0; i < dim_global_src[d]; i++) {
            ferr = fscanf(in,"%lf",&xg_src[i+offset]);
            if (ferr != 1) {
              printf("Error in ReadArraySerial(): unable to read data. ferr=%d\n", ferr);
              exit(1);
            }
          }
          offset += dim_global_src[d];
        }

        /* read solution */
        for (i = 0; i < nvars; i++) {
          int done = 0; _ArraySetValue_(index,ndims,0);
          while (!done) {
            int p; _ArrayIndex1D_(ndims,dim_global_src,index,0,p);
            ferr = fscanf(in,"%lf",&ug_src[p*nvars+i]);
            if (ferr != 1) {
              printf("Error in ReadArraySerial(): unable to read data. ferr=%d\n", ferr);
              exit(1);
            }
            _ArrayIncrementIndex_(ndims,dim_global_src,index,done);
          }
        }

        fclose(in);
      }
    } else if ((!strcmp(solver->ip_file_type,"bin")) || (!strcmp(solver->ip_file_type,"binary"))) {

      char filename[_MAX_STRING_SIZE_];
      strcpy(filename,fname_root);
      strcat(filename,".inp");
      FILE *in; in = fopen(filename,"rb");
      if (!in) *read_flag = 0;
      else {
        *read_flag = 1;
        printf("Reading array from binary file %s (Serial mode).\n",filename);
        size_t bytes;
        int size;
        /* allocate global solution array */
        size   = 1; for (d=0; d<ndims; d++) size *= dim_global_src[d]; size *= nvars;
        ug_src = (double*) calloc(size,sizeof(double));
        size   = 0; for (d=0; d<ndims; d++) size += dim_global_src[d];
        xg_src = (double*) calloc(size,sizeof(double));

        /* read grid */
        size = 0;
        for (d = 0; d < ndims; d++) size += dim_global_src[d];
        bytes = fread(xg_src, sizeof(double), size, in);
        if ((int)bytes != size) {
          fprintf(stderr,"Error in ReadArray(): Unable to read grid. Expected %d, Read %d.\n",
                  size, (int)bytes);
        }

        /* read solution */
        size = 1;
        for (d = 0; d < ndims; d++) size *= dim_global_src[d]; size *= nvars;
        bytes = fread(ug_src, sizeof(double), size, in);
        if ((int)bytes != size) {
          fprintf(stderr,"Error in ReadArray(): Unable to read solution. Expected %d, Read %d.\n",
                  size, (int)bytes);
        }

        fclose(in);
      }

    }

    double *ug_src_wg;
    long size = nvars;
    for (d=0; d<ndims; d++) size *= (dim_global_src[d]+2*ghosts); 
    ug_src_wg = (double*) calloc (size, sizeof(double));
    ArrayCopynD(  ndims, 
                  ug_src, 
                  ug_src_wg, 
                  dim_global_src, 
                  0, 
                  ghosts, 
                  index, 
                  nvars );
    fillGhostCells( dim_global_src, 
                    ghosts, 
                    ug_src_wg, 
                    nvars, 
                    ndims, 
                    solver->isPeriodic);
    free(ug_src);

    /* interpolate from the data read in to a global array 
     * with specified dimensions */
    int ierr = InterpolateGlobalnDVar(  dim_global,
                                        &ug,
                                        dim_global_src,
                                        ug_src_wg,
                                        nvars,
                                        ghosts,
                                        ndims,
                                        solver->isPeriodic );
    if (ierr) {
      fprintf(stderr, "Error in ReadArraywInterpSerial()\n");
      fprintf(stderr, "  InterpolateGlobalnDVar() returned with error!\n");
      return ierr;
    }

    if (x) {
      fprintf(stderr,"Error in ReadArraywInterpSerial()\n");
      fprintf(stderr,"  Not yet implemented for x\n");
      //InterpolateGlobal1DVar(xg_src, dim_global_src, xg, dim_global);
    }
    if (xg_src) free(xg_src);

    if (ug == NULL) {
      fprintf(stderr,"Error in ReadArraywInterp() - ug is NULL!\n");
      return 1;
    }
    if ((xg == NULL) && (x != NULL)) {
      fprintf(stderr,"Error in ReadArraywInterp() - xg is NULL!\n");
      return 1;
    }
  }

  /* Broadcast read_flag to all processes */
  IERR MPIBroadcast_integer(read_flag,1,0,&mpi->world); CHECKERR(ierr);

  if (*read_flag) {

    /* partition global array to all processes */
    IERR MPIPartitionArraynDwGhosts(  ndims,
                                      mpi,
                                      (mpi->rank?NULL:ug),
                                      u,dim_global,
                                      dim_local,
                                      ghosts,
                                      nvars ); CHECKERR(ierr);

//    if (x) {
//      /* partition x vector across the processes */
//      int offset_global = 0, offset_local = 0;
//      for (d=0; d<ndims; d++) {
//        IERR MPIPartitionArray1D(mpi,(mpi->rank?NULL:&xg[offset_global]),
//                                 &x[offset_local+ghosts],
//                                 mpi->is[d],mpi->ie[d],dim_local[d],0); CHECKERR(ierr);
//        offset_global += dim_global[d];
//        offset_local  += dim_local [d] + 2*ghosts;
//      }
//    }

    /* free global arrays */
    if (!mpi->rank) {
      free(ug);
      free(xg);
    }
  }

  return(0);
}