HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
Functions
FillGhostCells.c File Reference
#include <arrayfunctions.h>

Go to the source code of this file.

Functions

void fillGhostCells (const int *const a_dim, const int a_ngpt, double *const a_u, const int a_nvars, const int a_ndims, const int *const a_periodic)
 

Function Documentation

◆ fillGhostCells()

void fillGhostCells ( const int *const  a_dim,
const int  a_ngpt,
double *const  a_u,
const int  a_nvars,
const int  a_ndims,
const int *const  a_periodic 
)

Fill the ghost cells of a solution. Note that the solution array must be a global one (not one that is partitioned among MPI ranks). This is not a parallel operation (it will execute independently on multiple MPI ranks, if called by multiple processes).

For periodicity along any dimension, the ghost cells are filled appropriately from the other side of the domain. Otherwise, the interior data is extrapolated by a 4th order polynomial (assuming uniform grid spacing).

Parameters
a_dimgrid dimensions of solution
a_ngptnumber of ghost cells
a_usolution array
a_nvarsnumber of vector components of the solution
a_ndimsnumber of spatial dimensions
a_periodicperiodic or not along each dimension

Definition at line 17 of file FillGhostCells.c.

24 {
25  for (int d = 0; d < a_ndims; d++) {
26 
27  int bounds[a_ndims];
28  _ArrayCopy1D_(a_dim, bounds, a_ndims);
29  bounds[d] = a_ngpt;
30 
31  int index[a_ndims];
32  _ArraySetValue_(index, a_ndims, 0);
33 
34  if (a_periodic[d]) {
35 
36  /* periodic case */
37 
38  int done = 0;
39  while (!done) {
40 
41  {
42  /* low end - face = 1 */
43 
44  int p_gpt = 0,
45  p_int = 0;
46 
47  int index_gpt[a_ndims];
48  _ArrayCopy1D_(index, index_gpt, a_ndims);
49  index_gpt[d] -= a_ngpt;
50  _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
51 
52  int index_int[a_ndims];
53  _ArrayCopy1D_(index, index_int, a_ndims);
54  index_int[d] += (a_dim[d]-a_ngpt);
55  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int);
56 
57  _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
58  }
59 
60  {
61  /* high end - face = -1 */
62 
63  int p_gpt = 0,
64  p_int = 0;
65 
66  int index_gpt[a_ndims];
67  _ArrayCopy1D_(index, index_gpt, a_ndims);
68  index_gpt[d] += a_dim[d];
69  _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
70 
71  int index_int[a_ndims];
72  _ArrayCopy1D_(index, index_int, a_ndims);
73  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int);
74 
75  _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
76  }
77 
78  _ArrayIncrementIndex_(a_ndims, bounds, index, done);
79 
80  }
81 
82  } else {
83 
84  /* not periodic - extrapolate */
85 
86  int done = 0;
87  while (!done) {
88 
89  {
90  /* low end - face = 1 */
91 
92  int p_gpt = 0,
93  p_int_0 = 0,
94  p_int_1 = 0,
95  p_int_2 = 0,
96  p_int_3 = 0;
97 
98  int index_gpt[a_ndims];
99  _ArrayCopy1D_(index, index_gpt, a_ndims);
100  index_gpt[d] -= a_ngpt;
101  _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
102 
103  int index_int[a_ndims];
104  _ArrayCopy1D_(index, index_int, a_ndims);
105 
106  index_int[d] = 0;
107  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_0);
108  index_int[d]++;
109  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_1);
110  index_int[d]++;
111  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_2);
112  index_int[d]++;
113  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_3);
114 
115  double alpha = - (double) (a_ngpt - index[d]);
116  double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
117  double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
118  double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
119  double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;
120 
121  for (int v = 0; v < a_nvars; v++) {
122 
123  a_u[p_gpt*a_nvars+v] = c0 * a_u[p_int_0*a_nvars+v]
124  + c1 * a_u[p_int_1*a_nvars+v]
125  + c2 * a_u[p_int_2*a_nvars+v]
126  + c3 * a_u[p_int_3*a_nvars+v];
127 
128  }
129 
130  }
131 
132  {
133  /* high end - face = -1 */
134 
135  int p_gpt = 0,
136  p_int_0 = 0,
137  p_int_1 = 0,
138  p_int_2 = 0,
139  p_int_3 = 0;
140 
141  int index_gpt[a_ndims];
142  _ArrayCopy1D_(index, index_gpt, a_ndims);
143  index_gpt[d] += a_dim[d];
144  _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
145 
146  int index_int[a_ndims];
147  _ArrayCopy1D_(index, index_int, a_ndims);
148 
149  index_int[d] = a_dim[d]-1;
150  _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_0);
151  index_int[d]--;
152  _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_1);
153  index_int[d]--;
154  _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_2);
155  index_int[d]--;
156  _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_3);
157 
158  double alpha = - (double) (index[d]+1);
159  double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
160  double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
161  double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
162  double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;
163 
164  for (int v = 0; v < a_nvars; v++) {
165 
166  a_u[p_gpt*a_nvars+v] = c0 * a_u[p_int_0*a_nvars+v]
167  + c1 * a_u[p_int_1*a_nvars+v]
168  + c2 * a_u[p_int_2*a_nvars+v]
169  + c3 * a_u[p_int_3*a_nvars+v];
170 
171  }
172 
173  }
174 
175  _ArrayIncrementIndex_(a_ndims, bounds, index, done);
176 
177  }
178 
179  }
180 
181  }
182 
183  return;
184 }
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_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