HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
SparseGridsInterpolateGrid.cpp
Go to the documentation of this file.
1 
6 #include <cmath>
7 #include <basic.h>
8 #include <arrayfunctions.h>
9 #include <std_vec_ops.h>
11 
21  const SimulationObject* const a_src
22  )
23 {
24  /* get the source and destination grid dimensions */
25  GridDimensions dim_src, dim_dst;
26  StdVecOps::copyFrom(dim_src, a_src->solver.dim_global, m_ndims);
27  StdVecOps::copyFrom(dim_dst, a_dst->solver.dim_global, m_ndims);
28 
29  /* gather the source on rank 0 */
30  double* xg_src = NULL;
31  if (!m_rank) allocateGridArrays(dim_src, &xg_src);
32  {
33  int offset_global, offset_local;
34  offset_global = offset_local = 0;
35  for (int d=0; d<m_ndims; d++) {
36  MPIGatherArray1D( (void*) &(a_src->mpi),
37  (a_src->mpi.rank ? NULL : &xg_src[offset_global]),
38  &(a_src->solver.x[offset_local+a_src->solver.ghosts]),
39  a_src->mpi.is[d],
40  a_src->mpi.ie[d],
41  a_src->solver.dim_local[d],
42  0 );
43  offset_global += a_src->solver.dim_global[d];
44  offset_local += a_src->solver.dim_local [d] + 2*a_src->solver.ghosts;
45  }
46  }
47 
48  /* now do the interpolation, dimension-by-dimension */
49  double *xg_dst = NULL;
50  if (!m_rank) {
51 
52  GridDimensions dim_to(m_ndims,0);
53  GridDimensions dim_from(m_ndims,0);
54 
55  double* x_from;
56  double* x_to;
57 
58  dim_to = dim_src;
59  x_to = xg_src;
60  x_from = NULL;
61 
62  for (int dir = 0; dir < m_ndims; dir++) {
63 
64  dim_from = dim_to;
65  dim_to[dir] = dim_dst[dir];
66 
67  if (dim_from[dir] == dim_to[dir]) continue;
68 
69  double fac = (dim_to[dir] > dim_from[dir] ?
70  (double)dim_to[dir]/(double)dim_from[dir]
71  : (double)dim_from[dir]/(double)dim_to[dir] );
72  if (!isPowerOfTwo((int)fac)) {
73  fprintf(stderr,"Error in SparseGridsSimulation::interpolate() - \n");
74  fprintf(stderr," refinement/coarsening factor not a power of 2!\n");
75  exit(1);
76  }
77 
78  if (x_from != NULL) free(x_from);
79  x_from = x_to;
80 
81  allocateGridArrays(dim_to, &x_to);
82  if (dim_to[dir] < dim_from[dir]) {
83  coarsenGrid1D(dim_from, dim_to, x_from, x_to, dir);
84  } else {
85  refineGrid1D(dim_from, dim_to, x_from, x_to, dir);
86  }
87 
88  }
89 
90  /* dim_to should be equal to dim_dst now */
91  for (int d = 0; d < m_ndims; d++) {
92  if (dim_to[d] != dim_dst[d]) {
93  fprintf(stderr,"Error in SparseGridsSimulation::interpolate() - \n");
94  fprintf(stderr," dim_to[%d] != dim_dst[%d]!\n", d, d);
95  exit(1);
96  }
97  }
98 
99  if (x_from != NULL) free(x_from);
100  xg_dst = x_to;
101 
102  }
103 
104  /* partition the destination */
105  int offset_global = 0, offset_local = 0;
106  for (int d=0; d<m_ndims; d++) {
107  MPIPartitionArray1D( (void*) &(a_dst->mpi),
108  (a_dst->mpi.rank ? NULL : &xg_dst[offset_global]),
109  &(a_dst->solver.x[offset_local+a_dst->solver.ghosts]),
110  a_dst->mpi.is[d],
111  a_dst->mpi.ie[d],
112  a_dst->solver.dim_local[d],
113  0);
114  offset_global += a_dst->solver.dim_global[d];
115  offset_local += a_dst->solver.dim_local [d] + 2*a_dst->solver.ghosts;
116  }
117 
118  if (!m_rank) {
119  free(xg_dst);
120  }
121 
122  /* exchange MPI-boundary values of x between processors */
123  {
124  int offset = 0;
125  for (int d = 0; d < m_ndims; d++) {
126  MPIExchangeBoundaries1D( (void*) &(a_dst->mpi),
127  &(a_dst->solver.x[offset]),
128  a_dst->solver.dim_local[d],
129  a_dst->solver.ghosts,
130  d,
131  m_ndims);
132  offset += (a_dst->solver.dim_local[d] + 2*a_dst->solver.ghosts);
133  }
134  }
135  /* fill in ghost values of x at physical boundaries by extrapolation */
136  {
137  int offset = 0;
138  for (int d = 0; d < m_ndims; d++) {
139  double* X = &(a_dst->solver.x[offset]);
140  int* dim = a_dst->solver.dim_local;
141  int ghosts = a_dst->solver.ghosts;
142  if (a_dst->mpi.ip[d] == 0) {
143  /* fill left boundary along this dimension */
144  for (int i = 0; i < ghosts; i++) {
145  int delta = ghosts - i;
146  X[i] = X[ghosts] + ((double) delta) * (X[ghosts]-X[ghosts+1]);
147  }
148  }
149  if (a_dst->mpi.ip[d] == a_dst->mpi.iproc[d]-1) {
150  /* fill right boundary along this dimension */
151  for (int i = dim[d]+ghosts; i < dim[d]+2*ghosts; i++) {
152  int delta = i - (dim[d]+ghosts-1);
153  X[i] = X[dim[d]+ghosts-1]
154  + ((double) delta) * (X[dim[d]+ghosts-1]-X[dim[d]+ghosts-2]);
155  }
156  }
157  offset += (dim[d] + 2*ghosts);
158  }
159  }
160  return;
161 }
162 
169  const GridDimensions& a_dim_dst,
170  const double* const a_x_src,
171  double* const a_x_dst,
172  const int a_dir
173  )
174 {
175  for (int d = 0; d < m_ndims; d++) {
176  if ((d != a_dir) && (a_dim_src[d] != a_dim_dst[d])) {
177  fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
178  fprintf(stderr, " a_dim_src[%d] != a_dim_dst[%d]\n", d, d);
179  exit(1);
180  }
181  }
182 
183  int n_src = a_dim_src[a_dir];
184  int n_dst = a_dim_dst[a_dir];
185  if (n_dst > n_src) {
186  fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
187  fprintf(stderr, " destination grid is finer than source grid along a_dir!\n");
188  exit(1);
189  }
190 
191  double fac = ((double) n_src) / ((double) n_dst);
192  int stride = (int) fac;
193  if (std::abs(((double)stride)-fac) > _MACHINE_ZERO_) {
194  fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
195  fprintf(stderr, " non-integer coarsening factor!\n");
196  exit(1);
197  }
198 
199  const double* x_src = a_x_src;
200  double* x_dst = a_x_dst;
201 
202  for (int d = 0; d < m_ndims; d++) {
203 
204  if (d == a_dir) {
205 
206  for (int i = 0; i < a_dim_dst[d]; i++) {
207  double avg = 0;
208  for (int j=stride*i; j<stride*(i+1); j++) {
209  if (j >= n_src) {
210  fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
211  fprintf(stderr, " j >= n_src!\n");
212  exit(1);
213  }
214  avg += x_src[j];
215  }
216  avg /= (double) stride;
217  x_dst[i] = avg;
218  }
219 
220  } else {
221 
222  _ArrayCopy1D_(x_src, x_dst, a_dim_dst[d]);
223 
224  }
225 
226  x_src += a_dim_src[d];
227  x_dst += a_dim_dst[d];
228  }
229 
230  return;
231 }
232 
239  const GridDimensions& a_dim_dst,
240  const double* const a_x_src,
241  double* const a_x_dst,
242  const int a_dir
243  )
244 {
245  for (int d = 0; d < m_ndims; d++) {
246  if ((d != a_dir) && (a_dim_src[d] != a_dim_dst[d])) {
247  fprintf(stderr, "Error in SparseGridsSimulation::coarsenGrid1D() -\n");
248  fprintf(stderr, " a_dim_src[%d] != a_dim_dst[%d]\n", d, d);
249  exit(1);
250  }
251  }
252 
253  int n_src = a_dim_src[a_dir];
254  int n_dst = a_dim_dst[a_dir];
255  if (n_dst < n_src) {
256  fprintf(stderr, "Error in SparseGridsSimulation::refineGrid1D() -\n");
257  fprintf(stderr, " destination grid is coarser than source grid along a_dir!\n");
258  exit(1);
259  }
260 
261  double fac = ((double) n_dst) / ((double) n_src);
262  int stride = (int) fac;
263  if (std::abs(((double)stride)-fac) > _MACHINE_ZERO_) {
264  fprintf(stderr, "Error in SparseGridsSimulation::refineGrid1D() -\n");
265  fprintf(stderr, " non-integer refinement factor!\n");
266  exit(1);
267  }
268 
269  int offset = 0;
270  for (int d = 0; d < a_dir; d++) offset += a_dim_src[d];
271  const double* x_src = a_x_src + offset;
272  double* x_dst = a_x_dst + offset;
273 
274  fprintf(stderr, "Error in SparseGridsSimulation::refineGrid1D() -\n");
275  fprintf(stderr, " NOT YET IMPLEMENTED! Why do you need this?\n");
276  exit(1);
277 
278  return;
279 }
280 
Contains some vector ops.
int * dim_global
Definition: hypar.h:33
void copyFrom(std::vector< int > &a_iv, const int *const a_iv_carr, int a_n)
Definition: std_vec_ops.h:36
void interpolateGrid(SimulationObject *const, const SimulationObject *const)
int MPIGatherArray1D(void *, double *, double *, int, int, int, int)
int MPIExchangeBoundaries1D(void *, double *, int, int, int, int)
int MPIPartitionArray1D(void *, double *, double *, int, int, int, int)
int * dim_local
Definition: hypar.h:37
int ghosts
Definition: hypar.h:52
#define _ArrayCopy1D_(x, y, size)
#define _MACHINE_ZERO_
Definition: basic.h:26
void refineGrid1D(const GridDimensions &, const GridDimensions &, const double *const, double *const, int)
void coarsenGrid1D(const GridDimensions &, const GridDimensions &, const double *const, double *const, int)
void allocateGridArrays(const GridDimensions &a_dim, double **const a_x, const int a_ngpt=0)
Structure defining a simulation.
Some basic definitions and macros.
Contains macros and function definitions for common array operations.
double * x
Definition: hypar.h:107
std::vector< int > GridDimensions