HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
LinearADRAdvectionField.c
Go to the documentation of this file.
1 
6 #include <stdlib.h>
7 #include <math.h>
8 #include <string.h>
9 #include <basic.h>
10 #include <common.h>
11 #include <arrayfunctions.h>
12 #include <io.h>
14 #include <hypar.h>
15 #include <mpivars.h>
16 
34  void *m,
35  int idx,
36  int nsims,
37  int *dim_adv
38  )
39 {
40  HyPar *solver = (HyPar*) s;
41  MPIVariables *mpi = (MPIVariables*) m;
42  LinearADR *param = (LinearADR*) solver->physics;
43  double *adv = param->a;
44 
45  /* sanity checks */
46  if (param->constant_advection == 1) {
47  fprintf(stderr,"Error in LinearADRAdvectionField(): param->constant_advection == 1!\n");
48  return(1);
49  }
50  if (!strcmp(param->adv_filename,"none")) {
51  fprintf(stderr,"Error in LinearADRAdvectionField(): param->adv_filename is \"none\"!\n");
52  return(1);
53  }
54 
55  int d, done, flag;
56  int *dim = solver->dim_local;
57  int ghosts = solver->ghosts;
58 
59  /* number of advection field components at a grid point is
60  * number of spatial dimensions times number of solution
61  * components */
62  int adv_nvar = solver->ndims * solver->nvars;
63  if (param->adv_arr_size != adv_nvar*solver->npoints_local_wghosts) {
64  fprintf(stderr,"Error in LinearADRAdvectionField(): Incorrect adv_arr_size\n");
65  return(1);
66  }
67 
68  char fname_root[_MAX_STRING_SIZE_];
69  strcpy(fname_root, param->adv_filename);
70  if (idx >= 0) {
71  if (nsims > 1) {
72  char index[_MAX_STRING_SIZE_];
73  GetStringFromInteger(idx, index, (int)log10(nsims)+1);
74  strcat(fname_root, "_");
75  strcat(fname_root, index);
76  }
77  }
78 
79  /* read spatially-varying advection field from provided file */
80  if (dim_adv == NULL) {
81  int ierr = ReadArray( solver->ndims,
82  adv_nvar,
83  solver->dim_global,
84  solver->dim_local,
85  solver->ghosts,
86  solver,
87  mpi,
88  NULL,
89  adv,
90  fname_root,
91  &flag);
92  if (ierr) {
93  fprintf(stderr,"Error in LinearADRAdvectionField()\n");
94  fprintf(stderr," ReadArray() returned error!\n");
95  return ierr;
96  }
97  } else {
98  int ierr = ReadArraywInterp(solver->ndims,
99  adv_nvar,
100  solver->dim_global,
101  solver->dim_local,
102  dim_adv,
103  solver->ghosts,
104  solver,
105  mpi,
106  NULL,
107  adv,
108  fname_root,
109  &flag);
110  if (ierr) {
111  fprintf(stderr,"Error in LinearADRAdvectionField()\n");
112  fprintf(stderr," ReadArraywInterp() returned error!\n");
113  return ierr;
114  }
115  }
116 
117  if (!flag) {
118  /* if advection file not available, set it to zero */
119  _ArraySetValue_(adv,solver->npoints_local_wghosts*adv_nvar,0.0);
120  }
121 
122  /* if parallel, exchange MPI-boundary ghost point data */
124  adv_nvar,
125  solver->dim_local,
126  solver->ghosts,
127  mpi,
128  adv );
129 
130  /* Along each dimension:
131  - If boundaries are periodic, then set ghost points accordingly if
132  number of MPI ranks along that dimension is 1 (if more than 1,
133  MPIExchangeBoundariesnD() called above has taken care of it).
134  - Else, extrapolate at the boundaries. */
135  int indexb[solver->ndims],
136  indexi[solver->ndims],
137  bounds[solver->ndims],
138  offset[solver->ndims];
139  for (d = 0; d < solver->ndims; d++) {
140  if (solver->isPeriodic[d] && (mpi->iproc[d] == 1)) {
141  _ArrayCopy1D_(dim,bounds,solver->ndims); bounds[d] = ghosts;
142  /* left boundary */
143  done = 0; _ArraySetValue_(indexb,solver->ndims,0);
144  _ArraySetValue_(offset,solver->ndims,0); offset[d] = -ghosts;
145  while (!done) {
146  _ArrayCopy1D_(indexb,indexi,solver->ndims); indexi[d] = indexb[d] + dim[d] - ghosts;
147  int p1; _ArrayIndex1DWO_(solver->ndims,dim,indexb,offset,ghosts,p1);
148  int p2; _ArrayIndex1D_ (solver->ndims,dim,indexi,ghosts,p2);
149  _ArrayCopy1D_((adv+p2*adv_nvar), (adv+p1*adv_nvar), adv_nvar);
150  _ArrayIncrementIndex_(solver->ndims,bounds,indexb,done);
151  }
152  /* right boundary */
153  done = 0; _ArraySetValue_(indexb,solver->ndims,0);
154  _ArraySetValue_(offset,solver->ndims,0); offset[d] = dim[d];
155  while (!done) {
156  _ArrayCopy1D_(indexb,indexi,solver->ndims);
157  int p1; _ArrayIndex1DWO_(solver->ndims,dim,indexb,offset,ghosts,p1);
158  int p2; _ArrayIndex1D_ (solver->ndims,dim,indexi,ghosts,p2);
159  _ArrayCopy1D_((adv+p2*adv_nvar), (adv+p1*adv_nvar), adv_nvar);
160  _ArrayIncrementIndex_(solver->ndims,bounds,indexb,done);
161  }
162  } else {
163  /* left boundary */
164  if (!mpi->ip[d]) {
165  _ArrayCopy1D_(dim,bounds,solver->ndims); bounds[d] = ghosts;
166  _ArraySetValue_(offset,solver->ndims,0); offset[d] = -ghosts;
167  done = 0; _ArraySetValue_(indexb,solver->ndims,0);
168  while (!done) {
169  _ArrayCopy1D_(indexb,indexi,solver->ndims); indexi[d] = ghosts-1-indexb[d];
170  int p1; _ArrayIndex1DWO_(solver->ndims,dim,indexb,offset,ghosts,p1);
171  int p2; _ArrayIndex1D_ (solver->ndims,dim,indexi,ghosts,p2);
172  _ArrayCopy1D_((adv+p2*adv_nvar), (adv+p1*adv_nvar), adv_nvar);
173  _ArrayIncrementIndex_(solver->ndims,bounds,indexb,done);
174  }
175  }
176  /* right boundary */
177  if (mpi->ip[d] == mpi->iproc[d]-1) {
178  _ArrayCopy1D_(dim,bounds,solver->ndims); bounds[d] = ghosts;
179  _ArraySetValue_(offset,solver->ndims,0); offset[d] = dim[d];
180  done = 0; _ArraySetValue_(indexb,solver->ndims,0);
181  while (!done) {
182  _ArrayCopy1D_(indexb,indexi,solver->ndims); indexi[d] = dim[d]-1-indexb[d];
183  int p1; _ArrayIndex1DWO_(solver->ndims,dim,indexb,offset,ghosts,p1);
184  int p2; _ArrayIndex1D_ (solver->ndims,dim,indexi,ghosts,p2);
185  _ArrayCopy1D_((adv+p2*adv_nvar), (adv+p1*adv_nvar), adv_nvar);
186  _ArrayIncrementIndex_(solver->ndims,bounds,indexb,done);
187  }
188  }
189  }
190  }
191 
192  return(0);
193 }
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
char adv_filename[_MAX_STRING_SIZE_]
Definition: linearadr.h:43
int * dim_global
Definition: hypar.h:33
#define _ArrayIncrementIndex_(N, imax, i, done)
void GetStringFromInteger(int, char *, int)
int adv_arr_size
Definition: linearadr.h:47
void * physics
Definition: hypar.h:266
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
int * dim_local
Definition: hypar.h:37
MPI related function definitions.
double * a
Definition: linearadr.h:50
Function declarations for file I/O functions.
Linear Advection-Diffusion-Reaction model.
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
int ReadArray(int, int, int *, int *, int, void *, void *, double *, double *, char *, int *)
Definition: ReadArray.c:25
#define _MAX_STRING_SIZE_
Definition: basic.h:14
int * isPeriodic
Definition: hypar.h:162
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
int LinearADRAdvectionField(void *s, void *m, int idx, int nsims, int *dim_adv)
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
int ReadArraywInterp(int, int, int *, int *, int *, int, void *, void *, double *, double *, char *, int *)
Contains structure definition for hypar.
int constant_advection
Definition: linearadr.h:40
Some basic definitions and macros.
int ndims
Definition: hypar.h:26
Contains macros and function definitions for common array operations.
#define IERR
Definition: basic.h:16
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure of MPI-related variables.
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Some common functions used here and there.