HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ReadArraywInterp.c
Go to the documentation of this file.
1 
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <string.h>
9 #include <basic.h>
10 #include <arrayfunctions.h>
11 #include <mathfunctions.h>
12 #include <mpivars.h>
13 #include <hypar.h>
14 
15 static int ReadArraywInterpSerial(int,int,int*,int*,int*,int,void*,void*,double*,double*,char*,int*);
16 
28 int ReadArraywInterp( int ndims,
29  int nvars,
30  int *dim_global,
31  int *dim_local,
32  int *dim_global_src,
33  int ghosts,
34  void *s,
35  void *m,
36  double *x,
37  double *u,
38  char *fname_root,
39  int *read_flag
40  )
41 {
42  HyPar *solver = (HyPar*) s;
43  MPIVariables *mpi = (MPIVariables*) m;
45 
46  int retval = ReadArraywInterpSerial( ndims,
47  nvars,
48  dim_global,
49  dim_local,
50  dim_global_src,
51  ghosts,
52  s,
53  m,
54  x,
55  u,
56  fname_root,
57  read_flag );
58  if (retval) return retval;
59 
60  if (x) {
61  int offset, d;
62  /* exchange MPI-boundary values of x between processors */
63  offset = 0;
64  for (d = 0; d < ndims; d++) {
65  IERR MPIExchangeBoundaries1D(mpi,&x[offset],dim_local[d],
66  ghosts,d,ndims); CHECKERR(ierr);
67  offset += (dim_local [d] + 2*ghosts);
68  }
69  /* fill in ghost values of x at physical boundaries by extrapolation */
70  offset = 0;
71  for (d = 0; d < ndims; d++) {
72  double *X = &x[offset];
73  int *dim = dim_local, i;
74  if (mpi->ip[d] == 0) {
75  /* fill left boundary along this dimension */
76  for (i = 0; i < ghosts; i++) {
77  int delta = ghosts - i;
78  X[i] = X[ghosts] + ((double) delta) * (X[ghosts]-X[ghosts+1]);
79  }
80  }
81  if (mpi->ip[d] == mpi->iproc[d]-1) {
82  /* fill right boundary along this dimension */
83  for (i = dim[d]+ghosts; i < dim[d]+2*ghosts; i++) {
84  int delta = i - (dim[d]+ghosts-1);
85  X[i] = X[dim[d]+ghosts-1]
86  + ((double) delta) * (X[dim[d]+ghosts-1]-X[dim[d]+ghosts-2]);
87  }
88  }
89  offset += (dim[d] + 2*ghosts);
90  }
91  }
92 
93  return 0;
94 }
95 
159 int ReadArraywInterpSerial( int ndims,
160  int nvars,
161  int *dim_global,
162  int *dim_local,
163  int *dim_global_src,
164  int ghosts,
165  void *s,
166  void *m,
167  double *x,
168  double *u,
169  char *fname_root,
170  int *read_flag
171  )
172 {
173  HyPar *solver = (HyPar*) s;
174  MPIVariables *mpi = (MPIVariables*) m;
175  int i, d, ferr, index[ndims];
176  double *ug_src = NULL, *xg_src = NULL;
177  double *ug = NULL, *xg = NULL;
179 
180  *read_flag = 0;
181  /* Only root process reads from the file */
182  if (!mpi->rank) {
183 
184  /* read data from file - this data is of dimensions given by dim_global_sec */
185  if (!strcmp(solver->ip_file_type,"ascii")) {
186  char filename[_MAX_STRING_SIZE_];
187  strcpy(filename,fname_root);
188  strcat(filename,".inp");
189  FILE *in; in = fopen(filename,"r");
190  if (!in) *read_flag = 0;
191  else {
192  *read_flag = 1;
193  /* Reading from file */
194  printf("Reading array from ASCII file %s (Serial mode).\n",filename);
195  int size,offset;
196  /* allocate global solution array */
197  size = 1; for (d=0; d<ndims; d++) size *= dim_global_src[d]; size *= nvars;
198  ug_src = (double*) calloc(size,sizeof(double));
199  size = 0; for (d=0; d<ndims; d++) size += dim_global_src[d];
200  xg_src = (double*) calloc(size,sizeof(double));
201 
202  /* read grid */
203  offset = 0;
204  for (d = 0; d < ndims; d++) {
205  for (i = 0; i < dim_global_src[d]; i++) {
206  ferr = fscanf(in,"%lf",&xg_src[i+offset]);
207  if (ferr != 1) {
208  printf("Error in ReadArraySerial(): unable to read data. ferr=%d\n", ferr);
209  exit(1);
210  }
211  }
212  offset += dim_global_src[d];
213  }
214 
215  /* read solution */
216  for (i = 0; i < nvars; i++) {
217  int done = 0; _ArraySetValue_(index,ndims,0);
218  while (!done) {
219  int p; _ArrayIndex1D_(ndims,dim_global_src,index,0,p);
220  ferr = fscanf(in,"%lf",&ug_src[p*nvars+i]);
221  if (ferr != 1) {
222  printf("Error in ReadArraySerial(): unable to read data. ferr=%d\n", ferr);
223  exit(1);
224  }
225  _ArrayIncrementIndex_(ndims,dim_global_src,index,done);
226  }
227  }
228 
229  fclose(in);
230  }
231  } else if ((!strcmp(solver->ip_file_type,"bin")) || (!strcmp(solver->ip_file_type,"binary"))) {
232 
233  char filename[_MAX_STRING_SIZE_];
234  strcpy(filename,fname_root);
235  strcat(filename,".inp");
236  FILE *in; in = fopen(filename,"rb");
237  if (!in) *read_flag = 0;
238  else {
239  *read_flag = 1;
240  printf("Reading array from binary file %s (Serial mode).\n",filename);
241  size_t bytes;
242  int size;
243  /* allocate global solution array */
244  size = 1; for (d=0; d<ndims; d++) size *= dim_global_src[d]; size *= nvars;
245  ug_src = (double*) calloc(size,sizeof(double));
246  size = 0; for (d=0; d<ndims; d++) size += dim_global_src[d];
247  xg_src = (double*) calloc(size,sizeof(double));
248 
249  /* read grid */
250  size = 0;
251  for (d = 0; d < ndims; d++) size += dim_global_src[d];
252  bytes = fread(xg_src, sizeof(double), size, in);
253  if ((int)bytes != size) {
254  fprintf(stderr,"Error in ReadArray(): Unable to read grid. Expected %d, Read %d.\n",
255  size, (int)bytes);
256  }
257 
258  /* read solution */
259  size = 1;
260  for (d = 0; d < ndims; d++) size *= dim_global_src[d]; size *= nvars;
261  bytes = fread(ug_src, sizeof(double), size, in);
262  if ((int)bytes != size) {
263  fprintf(stderr,"Error in ReadArray(): Unable to read solution. Expected %d, Read %d.\n",
264  size, (int)bytes);
265  }
266 
267  fclose(in);
268  }
269 
270  }
271 
272  double *ug_src_wg;
273  long size = nvars;
274  for (d=0; d<ndims; d++) size *= (dim_global_src[d]+2*ghosts);
275  ug_src_wg = (double*) calloc (size, sizeof(double));
276  ArrayCopynD( ndims,
277  ug_src,
278  ug_src_wg,
279  dim_global_src,
280  0,
281  ghosts,
282  index,
283  nvars );
284  fillGhostCells( dim_global_src,
285  ghosts,
286  ug_src_wg,
287  nvars,
288  ndims,
289  solver->isPeriodic);
290  free(ug_src);
291 
292  /* interpolate from the data read in to a global array
293  * with specified dimensions */
294  int ierr = InterpolateGlobalnDVar( dim_global,
295  &ug,
296  dim_global_src,
297  ug_src_wg,
298  nvars,
299  ghosts,
300  ndims,
301  solver->isPeriodic );
302  if (ierr) {
303  fprintf(stderr, "Error in ReadArraywInterpSerial()\n");
304  fprintf(stderr, " InterpolateGlobalnDVar() returned with error!\n");
305  return ierr;
306  }
307 
308  if (x) {
309  fprintf(stderr,"Error in ReadArraywInterpSerial()\n");
310  fprintf(stderr," Not yet implemented for x\n");
311  //InterpolateGlobal1DVar(xg_src, dim_global_src, xg, dim_global);
312  }
313  if (xg_src) free(xg_src);
314 
315  if (ug == NULL) {
316  fprintf(stderr,"Error in ReadArraywInterp() - ug is NULL!\n");
317  return 1;
318  }
319  if ((xg == NULL) && (x != NULL)) {
320  fprintf(stderr,"Error in ReadArraywInterp() - xg is NULL!\n");
321  return 1;
322  }
323  }
324 
325  /* Broadcast read_flag to all processes */
326  IERR MPIBroadcast_integer(read_flag,1,0,&mpi->world); CHECKERR(ierr);
327 
328  if (*read_flag) {
329 
330  /* partition global array to all processes */
332  mpi,
333  (mpi->rank?NULL:ug),
334  u,dim_global,
335  dim_local,
336  ghosts,
337  nvars ); CHECKERR(ierr);
338 
339 // if (x) {
340 // /* partition x vector across the processes */
341 // int offset_global = 0, offset_local = 0;
342 // for (d=0; d<ndims; d++) {
343 // IERR MPIPartitionArray1D(mpi,(mpi->rank?NULL:&xg[offset_global]),
344 // &x[offset_local+ghosts],
345 // mpi->is[d],mpi->ie[d],dim_local[d],0); CHECKERR(ierr);
346 // offset_global += dim_global[d];
347 // offset_local += dim_local [d] + 2*ghosts;
348 // }
349 // }
350 
351  /* free global arrays */
352  if (!mpi->rank) {
353  free(ug);
354  free(xg);
355  }
356  }
357 
358  return(0);
359 }
360 
#define _ArraySetValue_(x, size, value)
static int ReadArraywInterpSerial(int, int, int *, int *, int *, int, void *, void *, double *, double *, char *, int *)
int MPIPartitionArraynDwGhosts(int, void *, double *, double *, int *, int *, int, int)
#define _ArrayIncrementIndex_(N, imax, i, done)
int MPIExchangeBoundaries1D(void *, double *, int, int, int, int)
INLINE int ArrayCopynD(int, const double *, double *, int *, int, int, int *, int)
MPI related function definitions.
void fillGhostCells(const int *const, const int, double *const, const int, const int, const int *const)
#define _ArrayIndex1D_(N, imax, i, ghost, index)
MPI_Comm world
#define _MAX_STRING_SIZE_
Definition: basic.h:14
int * isPeriodic
Definition: hypar.h:162
int InterpolateGlobalnDVar(const int *const, double **const, const int *const, double *const, const int, const int, const int, const int *const)
Contains function definitions for common mathematical functions.
char ip_file_type[_MAX_STRING_SIZE_]
Definition: hypar.h:180
int ReadArraywInterp(int, int, int *, int *, int *, int, void *, void *, double *, double *, char *, int *)
#define CHECKERR(ierr)
Definition: basic.h:18
Contains structure definition for hypar.
Some basic definitions and macros.
Contains macros and function definitions for common array operations.
#define IERR
Definition: basic.h:16
Structure of MPI-related variables.
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _DECLARE_IERR_
Definition: basic.h:17
int MPIBroadcast_integer(int *, int, int, void *)
Definition: MPIBroadcast.c:23