HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
SparseGridsFillGhostCells.cpp
Go to the documentation of this file.
1 
6 #include <arrayfunctions.h>
7 #include <sparse_grids_simulation.h>
8 
18 void SparseGridsSimulation::fillGhostCells( const GridDimensions& a_dim,
19  const int a_ngpt,
20  double* const a_u,
21  const int a_nvars
22  )
23 {
24  if (m_is_periodic.size() != m_ndims) {
25  fprintf(stderr, "Error in SparseGridsSimulation::fillGhostCells() -\n");
26  fprintf(stderr, "m_is_periodic.size() != m_ndims \n");
27  exit(1);
28  }
29 
30  for (int d = 0; d < m_ndims; d++) {
31 
32  int bounds[m_ndims];
33  _ArrayCopy1D_(a_dim, bounds, m_ndims);
34  bounds[d] = a_ngpt;
35 
36  int index[m_ndims];
37  _ArraySetValue_(index, m_ndims, 0);
38 
39  if (m_is_periodic[d]) {
40 
41  /* periodic case */
42 
43  int done = 0;
44  while (!done) {
45 
46  {
47  /* low end - face = 1 */
48 
49  int p_gpt = 0,
50  p_int = 0;
51 
52  int index_gpt[m_ndims];
53  _ArrayCopy1D_(index, index_gpt, m_ndims);
54  index_gpt[d] -= a_ngpt;
55  _ArrayIndex1D_(m_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
56 
57  int index_int[m_ndims];
58  _ArrayCopy1D_(index, index_int, m_ndims);
59  index_int[d] += (a_dim[d]-a_ngpt);
60  _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int);
61 
62  _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
63  }
64 
65  {
66  /* high end - face = -1 */
67 
68  int p_gpt = 0,
69  p_int = 0;
70 
71  int index_gpt[m_ndims];
72  _ArrayCopy1D_(index, index_gpt, m_ndims);
73  index_gpt[d] += a_dim[d];
74  _ArrayIndex1D_(m_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
75 
76  int index_int[m_ndims];
77  _ArrayCopy1D_(index, index_int, m_ndims);
78  _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int);
79 
80  _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
81  }
82 
83  _ArrayIncrementIndex_(m_ndims, bounds, index, done);
84 
85  }
86 
87  } else {
88 
89  /* not periodic - extrapolate */
90 
91  int done = 0;
92  while (!done) {
93 
94  {
95  /* low end - face = 1 */
96 
97  int p_gpt = 0,
98  p_int_0 = 0,
99  p_int_1 = 0,
100  p_int_2 = 0,
101  p_int_3 = 0;
102 
103  int index_gpt[m_ndims];
104  _ArrayCopy1D_(index, index_gpt, m_ndims);
105  index_gpt[d] -= a_ngpt;
106  _ArrayIndex1D_(m_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
107 
108  int index_int[m_ndims];
109  _ArrayCopy1D_(index, index_int, m_ndims);
110 
111  index_int[d] = 0;
112  _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int_0);
113  index_int[d]++;
114  _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int_1);
115  index_int[d]++;
116  _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int_2);
117  index_int[d]++;
118  _ArrayIndex1D_(m_ndims, a_dim, index_int, a_ngpt, p_int_3);
119 
120  double alpha = - (double) (a_ngpt - index[d]);
121  double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
122  double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
123  double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
124  double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;
125 
126  for (int v = 0; v < a_nvars; v++) {
127 
128  a_u[p_gpt*a_nvars+v] = c0 * a_u[p_int_0*a_nvars+v]
129  + c1 * a_u[p_int_1*a_nvars+v]
130  + c2 * a_u[p_int_2*a_nvars+v]
131  + c3 * a_u[p_int_3*a_nvars+v];
132 
133  }
134 
135  }
136 
137  {
138  /* high end - face = -1 */
139 
140  int p_gpt = 0,
141  p_int_0 = 0,
142  p_int_1 = 0,
143  p_int_2 = 0,
144  p_int_3 = 0;
145 
146  int index_gpt[m_ndims];
147  _ArrayCopy1D_(index, index_gpt, m_ndims);
148  index_gpt[d] += a_dim[d];
149  _ArrayIndex1D_(m_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
150 
151  int index_int[m_ndims];
152  _ArrayCopy1D_(index, index_int, m_ndims);
153 
154  index_int[d] = a_dim[d]-1;
155  _ArrayIndex1D_(m_ndims, a_dim, index, a_ngpt, p_int_0);
156  index_int[d]--;
157  _ArrayIndex1D_(m_ndims, a_dim, index, a_ngpt, p_int_1);
158  index_int[d]--;
159  _ArrayIndex1D_(m_ndims, a_dim, index, a_ngpt, p_int_2);
160  index_int[d]--;
161  _ArrayIndex1D_(m_ndims, a_dim, index, a_ngpt, p_int_3);
162 
163  double alpha = - (double) (index[d]+1);
164  double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
165  double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
166  double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
167  double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;
168 
169  for (int v = 0; v < a_nvars; v++) {
170 
171  a_u[p_gpt*a_nvars+v] = c0 * a_u[p_int_0*a_nvars+v]
172  + c1 * a_u[p_int_1*a_nvars+v]
173  + c2 * a_u[p_int_2*a_nvars+v]
174  + c3 * a_u[p_int_3*a_nvars+v];
175 
176  }
177 
178  }
179 
180  _ArrayIncrementIndex_(m_ndims, bounds, index, done);
181 
182  }
183 
184  }
185 
186  }
187 
188  return;
189 }
GridDimensions
std::vector< int > GridDimensions
Definition: sparse_grids_simulation.h:20
SparseGridsSimulation::fillGhostCells
void fillGhostCells(const GridDimensions &, const int, double *const, const int)
Definition: SparseGridsFillGhostCells.cpp:18
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
SparseGridsSimulation::m_is_periodic
std::vector< bool > m_is_periodic
Definition: sparse_grids_simulation.h:222
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
arrayfunctions.h
Contains macros and function definitions for common array operations.