cuIBM
A GPU-based Immersed Boundary Method code
bodies.cu
Go to the documentation of this file.
1 
7 #include "bodies.h"
8 
9 
16 template <typename memoryType>
17 void bodies<memoryType>::initialise(parameterDB &db, domain &D)
18 {
19  std::vector<body> *B = db["flow"]["bodies"].get<std::vector<body> *>();
20 
21  // number of bodies in the flow
22  numBodies = B->size();
23 
24  // set the sizes of all the arrays
25  numPoints.resize(numBodies);
26  offsets.resize(numBodies);
27 
28  startI.resize(numBodies);
29  startJ.resize(numBodies);
30  numCellsX.resize(numBodies);
31  numCellsY.resize(numBodies);
32 
33  xmin.resize(numBodies);
34  xmax.resize(numBodies);
35  ymin.resize(numBodies);
36  ymax.resize(numBodies);
37 
38  forceX.resize(numBodies);
39  forceY.resize(numBodies);
40 
41  // calculate offsets, number of points in each body and the total number of points
42  totalPoints = 0;
43  for(int k=0; k<numBodies; k++)
44  {
45  offsets[k] = totalPoints;
46  numPoints[k] = (*B)[k].numPoints;
47  totalPoints += numPoints[k];
48  }
49 
50  // fill up coordinates of body points
51  X.resize(totalPoints);
52  Y.resize(totalPoints);
53  ds.resize(totalPoints);
54  ones.resize(totalPoints);
55  cusp::blas::fill(ones, 1.0);
56  for(int k=0; k<numBodies; k++)
57  {
58  for(int i=0; i<numPoints[k]; i++)
59  {
60  X[i+offsets[k]] = (*B)[k].X[i];
61  Y[i+offsets[k]] = (*B)[k].Y[i];
62  }
63  }
64  x.resize(totalPoints);
65  y.resize(totalPoints);
66  uB.resize(totalPoints);
67  vB.resize(totalPoints);
68  I.resize(totalPoints);
69  J.resize(totalPoints);
70 
71  bodiesMove = false;
72  for(int k=0; k<numBodies; k++)
73  {
74  // assume a closed body (closed loop)
75  for(int i=offsets[k], j = offsets[k]+numPoints[k]-1; i<offsets[k]+numPoints[k];)
76  {
77  // calculate the lengths of the boundary segments
78  ds[i] = sqrt( (X[i]-X[j])*(X[i]-X[j]) + (Y[i]-Y[j])*(Y[i]-Y[j]) );
79 
80  // j takes the value of i, then i is incremented
81  j = i++;
82  }
83  // if the body is moving, set bodiesMove to true
84  bodiesMove = bodiesMove || (*B)[k].moving[0] || (*B)[k].moving[1];
85  }
86  // set initial position of the body
87  update(db, D, 0.0);
88 
89  if(numBodies)
90  {
91  calculateCellIndices(D);
92  calculateBoundingBoxes(db, D);
93  }
94 } // initialise
95 
96 
107 template <typename memoryType>
108 void bodies<memoryType>::calculateCellIndices(domain &D)
109 {
110  int i=0, j=0;
111 
112  // find the cell for the zeroth point
113  while(D.x[i+1] < x[0])
114  i++;
115  while(D.y[j+1] < y[0])
116  j++;
117  I[0] = i;
118  J[0] = j;
119 
120  for(int k=1; k<totalPoints; k++)
121  {
122  // if the next boundary point is to the left of the current boundary point
123  if(x[k] < x[k-1])
124  {
125  while(D.x[i] > x[k])
126  i--;
127  }
128  // if the next boundary point is to the right of the current boundary point
129  else
130  {
131  while(D.x[i+1] < x[k])
132  i++;
133  }
134  // if the next boundary point is below the current boundary point
135  if(y[k] < y[k-1])
136  {
137  while(D.y[j] > y[k])
138  j--;
139  }
140  // if the next boundary point is above the current boundary point
141  else
142  {
143  while(D.y[j+1] < y[k])
144  j++;
145  }
146  I[k] = i;
147  J[k] = j;
148  }
149 } // calculateCellIndices
150 
151 
163 template <typename memoryType>
164 void bodies<memoryType>::calculateBoundingBoxes(parameterDB &db, domain &D)
165 {
166  real scale = db["simulation"]["scaleCV"].get<real>(),
167  dx, dy;
168  int i, j;
169  for(int k=0; k<numBodies; k++)
170  {
171  xmin[k] = x[offsets[k]];
172  xmax[k] = xmin[k];
173  ymin[k] = y[offsets[k]];
174  ymax[k] = ymin[k];
175  for(int l=offsets[k]+1; l<offsets[k]+numPoints[k]; l++)
176  {
177  if(x[l] < xmin[k]) xmin[k] = x[l];
178  if(x[l] > xmax[k]) xmax[k] = x[l];
179  if(y[l] < ymin[k]) ymin[k] = y[l];
180  if(y[l] > ymax[k]) ymax[k] = y[l];
181  }
182  dx = xmax[k]-xmin[k];
183  dy = ymax[k]-ymin[k];
184  xmax[k] += 0.5*dx*(scale-1.0);
185  xmin[k] -= 0.5*dx*(scale-1.0);
186  ymax[k] += 0.5*dy*(scale-1.0);
187  ymin[k] -= 0.5*dy*(scale-1.0);
188 
189  i=0; j=0;
190  while(D.x[i+1] < xmin[k])
191  i++;
192  while(D.y[j+1] < ymin[k])
193  j++;
194  startI[k] = i;
195  startJ[k] = j;
196 
197  while(D.x[i] < xmax[k])
198  i++;
199  while(D.y[j] < ymax[k])
200  j++;
201  numCellsX[k] = i - startI[k];
202  numCellsY[k] = j - startJ[k];
203  }
204 } // calculateBoundingBoxes
205 
206 
222 template <typename memoryType>
223 void bodies<memoryType>::update(parameterDB &db, domain &D, real Time)
224 {
225  typedef typename cusp::array1d<real, memoryType> Array;
226  typedef typename Array::iterator Iterator;
227  typedef cusp::array1d_view<Iterator> View;
228 
229  // views of the vectors that store the coordinates and velocities of all the body points
230  View XView, YView, xView, yView, onesView, uBView, vBView;
231 
232  // body data
233  std::vector<body> *B = db["flow"]["bodies"].get<std::vector<body> *>();
234 
235  for(int l=0; l<numBodies; l++)
236  {
237  // update the location and velocity of the body
238  (*B)[l].update(Time);
239 
240  // create the views for the current body
241  if(l < numBodies-1)
242  {
243  XView = View(X.begin()+offsets[l], X.begin()+offsets[l+1]);
244  YView = View(Y.begin()+offsets[l], Y.begin()+offsets[l+1]);
245  xView = View(x.begin()+offsets[l], x.begin()+offsets[l+1]);
246  yView = View(y.begin()+offsets[l], y.begin()+offsets[l+1]);
247  onesView = View(ones.begin()+offsets[l], ones.begin()+offsets[l+1]);
248  uBView = View(uB.begin()+offsets[l], uB.begin()+offsets[l+1]);
249  vBView = View(vB.begin()+offsets[l], vB.begin()+offsets[l+1]);
250  }
251  else
252  {
253  XView = View(X.begin()+offsets[l], X.end());
254  YView = View(Y.begin()+offsets[l], Y.end());
255  xView = View(x.begin()+offsets[l], x.end());
256  yView = View(y.begin()+offsets[l], y.end());
257  onesView = View(ones.begin()+offsets[l], ones.end());
258  uBView = View(uB.begin()+offsets[l], uB.end());
259  vBView = View(vB.begin()+offsets[l], vB.end());
260  }
261 
262  // update postitions
263  // x-coordinates
264  cusp::blas::axpbypcz( onesView, XView, onesView, xView, (*B)[l].Xc[0], cos((*B)[l].Theta), -(*B)[l].X0[0]*cos((*B)[l].Theta) );
265  cusp::blas::axpbypcz( xView, YView, onesView, xView, 1.0, -sin((*B)[l].Theta), (*B)[l].X0[1]*sin((*B)[l].Theta) );
266  // y-coordinates
267  cusp::blas::axpbypcz( onesView, XView, onesView, yView, (*B)[l].Xc[1], sin((*B)[l].Theta), -(*B)[l].X0[0]*sin((*B)[l].Theta) );
268  cusp::blas::axpbypcz( yView, YView, onesView, yView, 1.0, cos((*B)[l].Theta), -(*B)[l].X0[1]*cos((*B)[l].Theta) );
269 
270  // update velocities
271  // x-velocities
272  cusp::blas::axpbypcz(onesView, yView, onesView, uBView, (*B)[l].vel[0], -(*B)[l].angVel, (*B)[l].angVel*(*B)[l].Xc[1]);
273  // y-velocities
274  cusp::blas::axpbypcz(onesView, xView, onesView, vBView, (*B)[l].vel[1], (*B)[l].angVel, -(*B)[l].angVel*(*B)[l].Xc[0]);
275  }
276 
277  if(numBodies)
278  calculateCellIndices(D);
279 } // update
280 
281 
288 template <>
289 void bodies<device_memory>::writeToFile(std::string &caseFolder, int timeStep)
290 {
291  vecH xHost = x,
292  yHost = y;
293  real *bx = thrust::raw_pointer_cast(&(xHost[0])),
294  *by = thrust::raw_pointer_cast(&(yHost[0]));
295  writeToFile(bx, by, caseFolder, timeStep);
296 } // writeToFile
297 
298 
307 template <typename memoryType>
308 void bodies<memoryType>::writeToFile(real *bx, real *by, std::string &caseFolder, int timeStep)
309 {
310  std::string path;
311  std::stringstream out;
312  out << caseFolder << '/' << std::setfill('0') << std::setw(7) << timeStep << "/bodies";
313  std::ofstream file(out.str().c_str());
314  file << '#' << std::setw(19) << "x-coordinate" << std::setw(20) << "y-coordinate" << std::endl;
315  for (int l=0; l < totalPoints; l++)
316  {
317  file << bx[l] << '\t' << by[l] << '\n';
318  }
319  file.close();
320 } // writeToFile
321 
322 
323 // specialization of the class bodies
324 template class bodies<device_memory>;
domain::x
vecH x
x-coordinates of the nodes
Definition: domain.h:22
real
double real
Is a float or a double depending on the machine precision.
Definition: types.h:116
bodies.h
Declaration of the class bodies.
bodies::calculateBoundingBoxes
void calculateBoundingBoxes(parameterDB &db, domain &D)
Calculates indices of the bounding box of each body in the flow.
Definition: bodies.cu:164
parameterDB
std::map< std::string, componentParameter > parameterDB
Map from a string to a componentParameter.
Definition: parameterDB.h:64
kernels::forceY
__global__ void forceY()
Kernel not usable.
Definition: calculateForce.cu:354
vecH
cusp::array1d< real, host_memory > vecH
Cusp 1D array stored in the host.
Definition: types.h:154
domain
Stores information about the computational grid.
Definition: domain.h:16
bodies
Contains information about bodies in the flow.
Definition: bodies.h:19
bodies::calculateCellIndices
void calculateCellIndices(domain &D)
Stores index of each cell that contains a boundary point.
Definition: bodies.cu:108
bodies::writeToFile
void writeToFile(std::string &caseFolder, int timeStep)
bodies::initialise
void initialise(parameterDB &db, domain &D)
Sets initial position and velocity of each body.
Definition: bodies.cu:17
bodies::update
void update(parameterDB &db, domain &D, real Time)
Updates position, velocity and neighbors of each body.
Definition: bodies.cu:223
kernels::forceX
__global__ void forceX(real *f, real *q, real *rn, int *tags, int nx, int ny, real *dx, real *dy, real dt, real alpha, real nu)
Kernel not usable.
Definition: calculateForce.cu:299
domain::y
vecH y
y-coordinates of the nodes
Definition: domain.h:22
Generated on Wed Jul 19 2017 12:27:42 for cuIBM by   doxygen 1.8.11