cuIBM
A GPU-based Immersed Boundary Method code
generateRN.cu
Go to the documentation of this file.
1 
7 #include "generateRN.h"
8 
9 
10 #define BSZ 16
11 
12 
17 namespace kernels
18 {
19 
37 __global__
38 void convectionTermU(real *rn, real *H, real *q,
39  int nx, int ny, real *dx, real *dy,
40  real dt, real gamma, real zeta, real alpha, real nu)
41 {
42  int bx = blockIdx.x,
43  by = blockIdx.y,
44  i = threadIdx.x,
45  j = threadIdx.y;
46 
47  // work out global index of first point in block
48  int I = (BSZ-2)*bx + i,
49  J = (BSZ-2)*by + j;
50 
51  if (I >= nx-1 || J >= ny) {
52  return;
53  }
54 
55  int N_u = (nx-1)*ny,
56  Gidx_x = J*(nx-1) + I,
57  Gidx_y = (J-1)*nx + I + N_u;
58 
59  real cTerm, dTerm, Hxn;
60 
61  __shared__ real u[BSZ][BSZ],
62  v[BSZ][BSZ],
63  Hx[BSZ][BSZ];
64 
65  __shared__ real Dx[BSZ][BSZ], Dy[BSZ][BSZ];
66 
67  Dy[j][i] = dy[J];
68  Dx[j][i] = dx[I];
69 
71  u[j][i] = q[Gidx_x]/Dy[j][i];
72  v[j][i] = q[Gidx_y]/Dx[j][i];
73  __syncthreads();
74 
76  int global_check = ( I==0 || I==(nx-2) || J==0 || J==(ny-1) ),
77  block_check = ( i==0 || i==(BSZ-1) || j==0 || j==(BSZ-1) );
78 
80  if( !(global_check || block_check) )
81  {
82  // copy Hx -> Hxn, Hy -> Hyn
83  Hxn = H[Gidx_x];
84 
85  // apply stencil
86  Hx[j][i] = - ( (u[j][i+1]+u[j][i])*(u[j][i+1]+u[j][i]) - (u[j][i]+u[j][i-1])*(u[j][i]+u[j][i-1]) )/( 2.0 * (Dx[j][i]+Dx[j][i+1]) ) \
87  - ( (u[j+1][i]+u[j][i])*(v[j+1][i]+v[j+1][i+1]) - (u[j][i]+u[j-1][i])*(v[j][i]+v[j][i+1]) )/( 4.0 * Dy[j][i] );
88 
89  H[Gidx_x] = Hx[j][i];
90 
91  // rN for u
92  cTerm = gamma*Hx[j][i] + zeta*Hxn;
93  //cTerm = Hx[j][i]; // 1st order Euler
94  dTerm = alpha*nu*2.0*( \
95  ( Dx[j][i]*u[j][i+1] - (Dx[j][i]+Dx[j][i+1])*u[j][i] + Dx[j][i+1]*u[j][i-1] )/( Dx[j][i]*Dx[j][i+1]*(Dx[j][i]+Dx[j][i+1]) ) \
96  + 4.0*( (Dy[j][i]+Dy[j-1][i])*u[j+1][i] - (Dy[j-1][i] + 2.0*Dy[j][i] + Dy[j+1][i])*u[j][i] + (Dy[j][i]+Dy[j+1][i])*u[j-1][i] ) \
97  /( (Dy[j][i]+Dy[j-1][i]) * (Dy[j-1][i] + 2.0*Dy[j][i] + Dy[j+1][i]) * (Dy[j][i]+Dy[j+1][i]) ) \
98  );
99  rn[Gidx_x] = (u[j][i]/dt + cTerm + dTerm) * 0.5*(Dx[j][i]+Dx[j][i+1]);
100  }
101 } // convectionTermU
102 
103 
121 __global__
122 void convectionTermV(real *rn, real *H, real *q,
123  int nx, int ny, real *dx, real *dy,
124  real dt, real gamma, real zeta, real alpha, real nu)
125 {
126  int bx = blockIdx.x,
127  by = blockIdx.y,
128  i = threadIdx.x,
129  j = threadIdx.y;
130 
131  // work out global index of first point in block
132  int I = (BSZ-2)*bx + i,
133  J = (BSZ-2)*by + j;
134 
135  if (I >= nx || J >= ny-1) {
136  return;
137  }
138 
139  int N_u = (nx-1)*ny,
140  Gidx_x = J*(nx-1) + I,
141  Gidx_y = J*nx + I + N_u;
142 
143  real cTerm, dTerm, Hyn;
144 
145  __shared__ real u[BSZ][BSZ],
146  v[BSZ][BSZ],
147  Hy[BSZ][BSZ];
148 
149  __shared__ real Dx[BSZ][BSZ], Dy[BSZ][BSZ];
150 
151  Dy[j][i] = dy[J];
152  Dx[j][i] = dx[I];
153 
155  u[j][i] = q[Gidx_x]/Dy[j][i];
156  v[j][i] = q[Gidx_y]/Dx[j][i];
157  __syncthreads();
158 
160  int global_check = ( I==0 || I==(nx-1) || J==0 || J==(ny-2) ),
161  block_check = ( i==0 || i==(BSZ-1) || j==0 || j==(BSZ-1) );
162 
164  if( !(global_check || block_check) )
165  {
166  // Y-component
167  // copy global data to the register (to store previous value)
168  Hyn = H[Gidx_y];
169 
170  // apply stencil
171  Hy[j][i] = - ( (u[j+1][i]+u[j][i])*(v[j][i]+v[j][i+1]) - (u[j+1][i-1]+u[j][i-1])*(v[j][i]+v[j][i-1]) )/(4.0 * Dx[j][i]) \
172  - ( (v[j+1][i]+v[j][i])*(v[j+1][i]+v[j][i]) - (v[j-1][i]+v[j][i])*(v[j-1][i]+v[j][i]) )/( 2.0 * (Dy[j+1][i]+Dy[j][i]) );
173 
174  // store newly calculated value in global memory (to be used in the next time step)
175  H[Gidx_y] = Hy[j][i];
176 
177  // rN for v
178  cTerm = gamma*Hy[j][i] + zeta*Hyn;
179  //cTerm = Hy[j][i]; // 1st order Euler
180  dTerm = alpha*nu*2.0*( \
181  4.0*( (Dx[j][i-1]+Dx[j][i])*v[j][i+1] - (Dx[j][i-1]+2.0*Dx[j][i]+Dx[j][i+1])*v[j][i] + (Dx[j][i]+Dx[j][i+1])*v[j][i-1] ) \
182  /( (Dx[j][i-1]+Dx[j][i]) * (Dx[j][i]+Dx[j][i+1]) * (Dx[j][i-1]+2.0*Dx[j][i]+Dx[j][i+1]) ) \
183  + ( Dy[j][i]*v[j+1][i] - (Dy[j][i]+Dy[j+1][i])*v[j][i] + Dy[j+1][i]*v[j-1][i] )/( Dy[j][i]*Dy[j+1][i]*(Dy[j][i]+Dy[j+1][i]) ) \
184  );
185  rn[Gidx_y] = (v[j][i]/dt + cTerm + dTerm) * 0.5*(Dy[j][i]+Dy[j+1][i]);
186  }
187 } // convectionTermV
188 
189 
210 __global__
211 void convectionTermVBottomTop(real *rn, real *H, real *q,
212  int nx, int ny, real *dx, real *dy,
213  real dt, real gamma, real zeta, real alpha, real nu,
214  real *bcBottom, real *bcTop)
215 {
216  int I = blockIdx.x*blockDim.x + threadIdx.x;
217 
218  if(I==0 || I >= nx-1) return;
220  int Iu, Iv = I + (nx-1)*ny;
221  real Hn, cTerm, dTerm;
223  Hn = H[Iv];
224  H[Iv] = -( \
225  ( 0.5*(q[Iv]/dx[I]+q[Iv+1]/dx[I+1]) ) * ( 0.5*(q[I]/dy[0]+q[I+(nx-1)]/dy[1]) ) \
226  - ( 0.5*(q[Iv-1]/dx[I-1]+q[Iv]/dx[I]) ) * ( 0.5*(q[I-1]/dy[0]+q[I-1+(nx-1)]/dy[1]) ) \
227  )/dx[I] \
228  -( \
229  (0.5*(q[Iv] + q[Iv+nx])/dx[I]) * (0.5*(q[Iv] + q[Iv+nx])/dx[I]) \
230  - (0.5*(bcBottom[I+nx-1] + q[Iv]/dx[I])) * (0.5*(bcBottom[I+nx-1] + q[Iv]/dx[I])) \
231  )/(0.5*(dy[0] + dy[1]));
232  cTerm = gamma*H[Iv] + zeta*Hn;
233  //cTerm = H[Iv]; /// 1st order Euler
235  dTerm = alpha*nu*2.0*( \
236  4.0 * ( (dx[I-1]+dx[I])*q[Iv+1]/dx[I+1] - (dx[I-1]+2.0*dx[I]+dx[I+1])*q[Iv]/dx[I] + (dx[I]+dx[I+1])*q[Iv-1]/dx[I-1] ) \
237  /( (dx[I-1]+dx[I]) * (dx[I]+dx[I+1]) * (dx[I-1]+2.0*dx[I]+dx[I+1]) ) \
238  + ( dy[0]*q[Iv+nx]/dx[I] - (dy[0]+dy[1])*q[Iv]/dx[I] + dy[1]*bcBottom[I+nx-1] )/( dy[0]*dy[1]*(dy[0]+dy[1]) ) \
239  );
241  rn[Iv] = ( q[Iv]/(dx[I]*dt) + cTerm + dTerm ) * 0.5*(dy[0] + dy[1]);
243  Iu = I + (ny-2)*(nx-1);
244  Iv = I + (nx-1)*ny + (ny-2)*nx;
246  Hn = H[Iv];
247  H[Iv] = -(
248  ( 0.5*(q[Iv]/dx[I] + q[Iv+1]/dx[I+1]) )*( 0.5*(q[Iu]/dy[ny-2] + q[Iu+(nx-1)]/dy[ny-1]) ) \
249  - ( 0.5*(q[Iv-1]/dx[I-1] + q[Iv]/dx[I]) )*( 0.5*(q[Iu-1]/dy[ny-2] + q[Iu-1+(nx-1)]/dy[ny-1]) )
250  )/(dx[I])
251  -(
252  ( 0.5*(q[Iv]/dx[I] + bcTop[I+nx-1]) )*( 0.5*(q[Iv]/dx[I] + bcTop[I+nx-1]) ) \
253  - ( 0.5*(q[Iv-nx] + q[Iv])/dx[I] )*( 0.5*(q[Iv-nx] + q[Iv])/dx[I] )
254  )/(0.5*(dy[ny-2] + dy[ny-1]));
255  cTerm = gamma*H[Iv] + zeta*Hn;
256  //cTerm = H[Iv]; /// 1st order Euler
258  dTerm = alpha*nu*2.0*( \
259  4.0*( (dx[I-1]+dx[I])*q[Iv+1]/dx[I+1] - (dx[I-1]+2.0*dx[I]+dx[I+1])*q[Iv]/dx[I] + (dx[I]+dx[I+1])*q[Iv-1]/dx[I-1] ) \
260  /( (dx[I-1]+dx[I]) * (dx[I]+dx[I+1]) * (dx[I-1]+2.0*dx[I]+dx[I+1]) ) \
261  + ( dy[ny-2]*bcTop[I+nx-1] - (dy[ny-1]+dy[ny-2])*q[Iv]/dx[I] + dy[ny-1]*q[Iv-nx]/dx[I] )/( dy[ny-2]*dy[ny-1]*(dy[ny-2]+dy[ny-1]) ) \
262  );
264  rn[Iv] = ( q[Iv]/(dx[I]*dt) + cTerm + dTerm ) * 0.5*(dy[ny-2] + dy[ny-1]);
265 } // convectionTermVBottomTop
266 
267 
290 __global__
291 void convectionTermUBottomTop(real *rn, real *H, real *q,
292  int nx, int ny, real *dx, real *dy,
293  real dt, real gamma, real zeta, real alpha, real nu,
294  real *bcBottom, real *bcTop, real *bcLeft, real *bcRight)
295 {
296  int I = blockIdx.x*blockDim.x + threadIdx.x;
297 
298  if(I >= nx-1) return;
300  int Iu,
301  Iv = I + (nx-1)*ny;
302  real u = q[I]/dy[0],
303  u0x, u1x,
304  ul, ur,
305  Hn, cTerm, dTerm;
307  if(I==0){
308  u0x = 0.5*(bcLeft[0] + u);
309  u1x = 0.5*(u + q[I+1]/dy[0]);
310  ul = bcLeft[0];
311  ur = q[I+1]/dy[0];
312  }
313  else if(I==nx-2){
314  u0x = 0.5*(q[I-1]/dy[0] + u);
315  u1x = 0.5*(u + bcRight[0]);
316  ul = q[I-1]/dy[0];
317  ur = bcRight[0];
318  }
319  else{
320  u0x = 0.5*(q[I-1]/dy[0] + u);
321  u1x = 0.5*(u + q[I+1]/dy[0]);
322  ul = q[I-1]/dy[0];
323  ur = q[I+1]/dy[0];
324  }
326  Hn = H[I];
327  H[I] = - ( u1x*u1x - u0x*u0x )/( 0.5*(dx[I]+dx[I+1]) ) \
328  - ( 0.5*(q[Iv]/dx[I] + q[Iv+1]/dx[I+1]) * 0.5*(u + q[I+(nx-1)]/dy[1]) - 0.5*(bcBottom[I+(nx-1)] + bcBottom[I+1+(nx-1)])*bcBottom[I] )/(dy[0]);
329  cTerm = gamma*H[I] + zeta*Hn;
330  //cTerm = H[I]; // 1st order Euler **************************** DID NOT CHANGE I TO Iu HERE
331  // Diffusion Term
332  dTerm = alpha*nu*2.0*( \
333  ( dx[I]*ur - (dx[I]+dx[I+1])*u + dx[I+1]*ul )/( dx[I] * (dx[I]+dx[I+1]) * dx[I+1] ) \
334  + 4.0*( dy[0]*q[I+(nx-1)]/dy[1] - (2.0*dy[0]+dy[1])*u + (dy[0]+dy[1])*bcBottom[I] ) \
335  /( dy[0] * (2.0*dy[0]+dy[1]) * (dy[0]+dy[1]) )
336  );
338  rn[I] = ( u/dt + cTerm + dTerm ) * 0.5*(dx[I] + dx[I+1]);
340  Iu = I + (ny-1)*(nx-1);
341  Iv = I + (ny-2)*nx + (nx-1)*ny;
342  u = q[Iu]/dy[ny-1];
344  if(I==0){
345  u0x = 0.5*(bcLeft[ny-1] + u);
346  u1x = 0.5*(u + q[Iu+1]/dy[ny-1]);
347  ul = bcLeft[ny-1];
348  ur = q[Iu+1]/dy[ny-1];
349  }
350  else if(I==nx-2){
351  u0x = 0.5*(q[Iu-1]/dy[ny-1] + u);
352  u1x = 0.5*(u + bcRight[ny-1]);
353  ul = q[Iu-1]/dy[ny-1];
354  ur = bcRight[ny-1];
355  }
356  else{
357  u0x = 0.5*(q[Iu-1]/dy[ny-1] + u);
358  u1x = 0.5*(u + q[Iu+1]/dy[ny-1]);
359  ul = q[Iu-1]/dy[ny-1];
360  ur = q[Iu+1]/dy[ny-1];
361  }
363  Hn = H[Iu];
364  H[Iu] = - ( u1x*u1x - u0x*u0x )/( 0.5*(dx[I]+dx[I+1]) ) \
365  - ( bcTop[I]*0.5*(bcTop[I+(nx-1)]+bcTop[I+1+(nx-1)]) - 0.5*(q[Iv]/dx[I] + q[Iv+1]/dx[I+1])*0.5*(u + q[Iu-(nx-1)]/dy[ny-2]) )/(dy[ny-1]);
366  cTerm = gamma*H[Iu] + zeta*Hn;
367  //cTerm = H[Iu]; /// 1st order Euler ************************* CHANGED I TO Iu HERE
368  // Diffusion Term
369  dTerm = alpha*nu*2.0*( \
370  ( dx[I]*ur - (dx[I]+dx[I+1])*u + dx[I+1]*ul )/( dx[I] * (dx[I]+dx[I+1]) * dx[I+1] ) \
371  + 4.0*( (dy[ny-2]+dy[ny-1])*bcTop[I] - (dy[ny-2]+2.0*dy[ny-1])*u + (dy[ny-1])*q[Iu-(nx-1)]/dy[ny-2] ) \
372  /( (dy[ny-2]+dy[ny-1])*(dy[ny-1])*(dy[ny-2]+2.0*dy[ny-1]) )
373  );
375  rn[Iu] = ( u/dt + cTerm + dTerm) * 0.5*(dx[I] + dx[I+1]);
376 } // convectionTermUBottomTop
377 
378 
399 __global__
400 void convectionTermULeftRight(real *rn, real *H, real *q,
401  int nx, int ny, real *dx, real *dy,
402  real dt, real gamma, real zeta, real alpha, real nu,
403  real *bcLeft, real *bcRight)
404 {
405  int I = blockIdx.x*blockDim.x + threadIdx.x;
406 
407  if(I==0 || I >= ny-1) return;
409  int Iu = I*(nx-1),
410  Iv = I*(nx) + (nx-1)*ny;
411  real Hn, cTerm, dTerm;
413  Hn = H[Iu];
414  H[Iu] = -( \
415  ( 0.5*(q[Iu] + q[Iu+1])/dy[I] ) * ( 0.5*(q[Iu] + q[Iu+1])/dy[I] ) \
416  - ( 0.5*(bcLeft[I] + q[Iu]/dy[I]) ) * ( 0.5*(bcLeft[I] + q[Iu]/dy[I]) ) \
417  )/(0.5*(dx[0]+dx[1])) \
418  -( \
419  ( 0.5*(q[Iv]/dx[0] + q[Iv+1]/dx[1]) ) * ( 0.5*(q[Iu]/dy[I] + q[Iu+(nx-1)]/dy[I+1]) ) \
420  - ( 0.5*(q[Iv-nx]/dx[0] + q[Iv+1-nx]/dx[1]) ) * ( 0.5*(q[Iu-(nx-1)]/dy[I-1] + q[Iu]/dy[I]) ) \
421  )/(dy[I]);
422  cTerm = gamma*H[Iu] + zeta*Hn;
423  //cTerm = H[Iu]; /// 1st order Euler
425  dTerm = alpha*nu*2.0*( \
426  ( dx[0]*q[Iu+1]/dy[I] - (dx[0]+dx[1])*q[Iu]/dy[I] + dx[1]*bcLeft[I] )/( dx[0]*dx[1]*(dx[0]+dx[1]) ) \
427  + 4.0 * ( (dy[I-1]+dy[I])*q[Iu+(nx-1)]/dy[I+1] - (dy[I-1]+2.0*dy[I]+dy[I+1])*q[Iu]/dy[I] + (dy[I]+dy[I+1])*q[Iu-(nx-1)]/dy[I-1] ) \
428  /( (dy[I-1]+dy[I]) * (dy[I]+dy[I+1]) * (dy[I-1]+2.0*dy[I]+dy[I+1]) ) \
429  );
431  rn[Iu] = ( q[Iu]/(dy[I]*dt) + cTerm + dTerm ) * 0.5*(dx[0] + dx[1]);
433  Iu = I*(nx-1) + (nx-2);
434  Iv = I*(nx) + (nx-1)*ny + (nx-2);
436  H[Iu] = -( \
437  ( 0.5*(q[Iu]/dy[I] + bcRight[I]) ) * ( 0.5*(q[Iu]/dy[I] + bcRight[I]) ) \
438  - ( 0.5*(q[Iu-1] + q[Iu])/dy[I] ) * ( 0.5*(q[Iu-1] + q[Iu])/dy[I] ) \
439  )/(0.5*(dx[nx-2]+dx[nx-1])) \
440  -( \
441  ( 0.5*(q[Iv]/dx[nx-2] + q[Iv+1]/dx[nx-1]) ) * ( 0.5*(q[Iu]/dy[I] + q[Iu+(nx-1)]/dy[I+1]) ) \
442  - ( 0.5*(q[Iv-nx]/dx[nx-2] + q[Iv+1-nx]/dx[nx-1]) ) * ( 0.5*(q[Iu-(nx-1)]/dy[I-1] + q[Iu]/dy[I]) ) \
443  )/(dy[I]);
444  cTerm = gamma*H[Iu] + zeta*Hn;
445  //cTerm = H[Iu]; /// 1st order Euler
447  dTerm = alpha*nu*2.0*( \
448  ( dx[nx-2]*bcRight[I] - (dx[nx-2]+dx[nx-1])*q[Iu]/dy[I] + dx[nx-1]*q[Iu-1]/dy[I] )/( dx[nx-2]*dx[nx-1]*(dx[nx-2]+dx[nx-1]) ) \
449  + 4.0 * ( (dy[I-1]+dy[I])*q[Iu+(nx-1)]/dy[I+1] - (dy[I-1]+2.0*dy[I]+dy[I+1])*q[Iu]/dy[I] + (dy[I]+dy[I+1])*q[Iu-(nx-1)]/dy[I-1] ) \
450  /( (dy[I-1]+dy[I]) * (dy[I]+dy[I+1]) * (dy[I-1]+2.0*dy[I]+dy[I+1]) ) \
451  );
453  rn[Iu] = ( q[Iu]/(dy[I]*dt) + cTerm + dTerm ) * 0.5*(dx[nx-2] + dx[nx-1]);
454 } // convectionTermULeftRight
455 
456 
479 __global__
480 void convectionTermVLeftRight(real *rn, real *H, real *q,
481  int nx, int ny, real *dx, real *dy,
482  real dt, real gamma, real zeta, real alpha, real nu,
483  real *bcBottom, real *bcTop, real *bcLeft, real *bcRight)
484 {
485  int I = blockIdx.x*blockDim.x + threadIdx.x;
486 
487  if(I > ny-2) return;
489  int Iu = I*(nx-1),
490  Iv = I*nx + (nx-1)*ny;
491  real vb, vt, v0y, v1y, v,
492  Hn, cTerm, dTerm;
493  v = q[Iv]/dx[0];
494  if(I==0){
495  v0y = 0.5*(bcBottom[nx-1] + v);
496  v1y = 0.5*(v + q[Iv+nx]/dx[0]);
497  vb = bcBottom[nx-1];
498  vt = q[Iv+nx]/dx[0];
499  }
500  else if(I==ny-2){
501  v0y = 0.5*(q[Iv-nx]/dx[0] + v);
502  v1y = 0.5*(q[Iv]/dx[0] + bcTop[nx-1]);
503  vb = q[Iv-nx]/dx[0];
504  vt = bcTop[nx-1];
505  }
506  else{
507  v0y = 0.5*(q[Iv-nx]/dx[0] + v);
508  v1y = 0.5*(v + q[Iv+nx]/dx[0]);
509  vb = q[Iv-nx]/dx[0];
510  vt = q[Iv+nx]/dx[0];
511  }
512  __syncthreads();
514  Hn = H[Iv];
515  H[Iv] = -( 0.5*(v + q[Iv+1]/dx[1])*0.5*(q[Iu]/dy[I] + q[Iu+(nx-1)]/dy[I+1]) - bcLeft[I+ny]*0.5*(bcLeft[I] + bcLeft[I+1]) )/dx[0] \
516  -( v1y*v1y - v0y*v0y )/(0.5*(dy[I] + dy[I+1]));
517  cTerm = gamma*H[Iv] + zeta*Hn;
518  //cTerm = H[Iv]; /// 1st order Euler
520  dTerm = alpha*nu*2.0*( \
521  4.0 * ( (dx[0])*q[Iv+1]/dx[1] - (2.0*dx[0]+dx[1])*v + (dx[0]+dx[1])*bcLeft[I+ny] ) \
522  /( (dx[0]) * (dx[0]+dx[1]) * (2.0*dx[0]+dx[1]) ) \
523  + ( dy[I]*vt - (dy[I]+dy[I+1])*v + dy[I+1]*vb )/( dy[I]*dy[I+1]*(dy[I]+dy[I+1]) ) \
524  );
526  rn[Iv] = ( v/dt + cTerm + dTerm ) * 0.5*(dy[I] + dy[I+1]);
528  Iu = I*(nx-1) + (nx-1);
529  Iv = I*nx + (nx-1)*ny + (nx-1);
530  v = q[Iv]/dx[nx-1];
531  if(I==0){
532  v0y = 0.5*(bcBottom[nx-1+(nx-1)] + v);
533  v1y = 0.5*(v + q[Iv+nx]/dx[nx-1]);
534  vb = bcBottom[nx-1+(nx-1)];
535  vt = q[Iv+nx]/dx[nx-1];
536  }
537  else if(I==ny-2){
538  v0y = 0.5*(q[Iv-nx]/dx[nx-1] + v);
539  v1y = 0.5*(q[Iv]/dx[nx-1] + bcTop[nx-1+(nx-1)]);
540  vb = q[Iv-nx]/dx[nx-1];
541  vt = bcTop[nx-1+(nx-1)];
542  }
543  else{
544  v0y = 0.5*(q[Iv-nx]/dx[nx-1] + v);
545  v1y = 0.5*(v + q[Iv+nx]/dx[nx-1]);
546  vb = q[Iv-nx]/dx[nx-1];
547  vt = q[Iv+nx]/dx[nx-1];
548  }
549  __syncthreads();
551  Hn = H[Iv];
552  H[Iv] = -( bcRight[I+ny]*0.5*(bcRight[I]+bcRight[I+1]) - 0.5*(q[Iv-1]/dx[nx-2] + v)*0.5*(q[Iu-1]/dy[I] + q[Iu-1+(nx-1)]/dy[I+1]) )/dx[nx-1] \
553  -( v1y*v1y - v0y*v0y )/(0.5*(dy[I] + dy[I+1]));
554  cTerm = gamma*H[Iv] + zeta*Hn;
555  //cTerm = H[Iv]; /// 1st order Euler
557  dTerm = alpha*nu*2.0*( \
558  4.0 * ( (dx[nx-1]+dx[nx-2])*bcRight[I+ny] - (2.0*dx[nx-1]+dx[nx-2])*v + (dx[nx-1])*q[Iv-1]/dx[nx-2] ) \
559  /( (dx[nx-1]) * (dx[nx-1]+dx[nx-2]) * (2.0*dx[nx-1]+dx[nx-2]) ) \
560  + ( dy[I]*vt - (dy[I]+dy[I+1])*v + dy[I+1]*vb )/( dy[I]*dy[I+1]*(dy[I]+dy[I+1]) ) \
561  );
563  rn[Iv] = ( v/dt + cTerm + dTerm ) * 0.5*(dy[I] + dy[I+1]);
564 } // convectionTermVLeftRight
565 
566 } // End of namespace kernels
generateRN.h
Declaration of the kernels to generate the explicit terms of the momentum equation.
real
double real
Is a float or a double depending on the machine precision.
Definition: types.h:116
kernels
Contains all the custom-written CUDA kernels.
kernels::convectionTermVLeftRight
__global__ void convectionTermVLeftRight(real *rn, real *H, real *q, int nx, int ny, real *dx, real *dy, real dt, real gamma, real zeta, real alpha, real nu, real *bcBottom, real *bcTop, real *bcLeft, real *bcRight)
Computes v-component of the vector resulting from the explicit terms of the discretized momentum equa...
Definition: generateRN.cu:480
kernels::convectionTermUBottomTop
__global__ void convectionTermUBottomTop(real *rn, real *H, real *q, int nx, int ny, real *dx, real *dy, real dt, real gamma, real zeta, real alpha, real nu, real *bcBottom, real *bcTop, real *bcLeft, real *bcRight)
Computes u-component of the vector resulting from the explicit terms of the discretized momentum equa...
Definition: generateRN.cu:291
BSZ
#define BSZ
Definition: generateRN.cu:10
kernels::convectionTermVBottomTop
__global__ void convectionTermVBottomTop(real *rn, real *H, real *q, int nx, int ny, real *dx, real *dy, real dt, real gamma, real zeta, real alpha, real nu, real *bcBottom, real *bcTop)
Computes v-component of the vector resulting from the explicit terms of the discretized momentum equa...
Definition: generateRN.cu:211
kernels::convectionTermU
__global__ void convectionTermU(real *rn, real *H, real *q, int nx, int ny, real *dx, real *dy, real dt, real gamma, real zeta, real alpha, real nu)
Computes u-component of the vector resulting from the explicit terms of the discretized momentum equa...
Definition: generateRN.cu:38
kernels::convectionTermV
__global__ void convectionTermV(real *rn, real *H, real *q, int nx, int ny, real *dx, real *dy, real dt, real gamma, real zeta, real alpha, real nu)
Computes v-component of the vector resulting from the explicit terms of the discretized momentum equa...
Definition: generateRN.cu:122
kernels::convectionTermULeftRight
__global__ void convectionTermULeftRight(real *rn, real *H, real *q, int nx, int ny, real *dx, real *dy, real dt, real gamma, real zeta, real alpha, real nu, real *bcLeft, real *bcRight)
Computes u-component of the vector resulting from the explicit terms of the discretized momentum equa...
Definition: generateRN.cu:400
Generated on Wed Jul 19 2017 12:27:42 for cuIBM by   doxygen 1.8.11