cuIBM
A GPU-based Immersed Boundary Method code
generateL.inl
Go to the documentation of this file.
1 
12 template <>
13 void DirectForcingSolver<device_memory>::generateL()
14 {
15 #if 1
16  int nx = domInfo->nx,
17  ny = domInfo->ny;
18 
19  real Cx0 = 0.0, Cx1 = 0.0, Cy0 = 0.0, Cy1 = 0.0,
20  scale = 0.0,
21  dx0 = 1.0, dx1 = 1.0, dy0 = 1.0, dy1 = 1.0;
22 
23  real *dx = thrust::raw_pointer_cast(&(domInfo->dx[0])),
24  *dy = thrust::raw_pointer_cast(&(domInfo->dy[0]));
25 
26  int numU = (nx-1)*ny;
27  int numUV = numU + nx*(ny-1);
28  int I,
29  num_elements = 0;
30 
31  cooH LHost(numUV, numUV, 5*numUV-4*(nx+ny)+4);
32 
33  real nu = (*paramDB)["flow"]["nu"].get<real>();
34  real dt = (*paramDB)["simulation"]["dt"].get<real>();
35  interpolationType interpType = (*paramDB)["simulation"]["interpolationType"].get<interpolationType>();
36 
37  bool notBdryNode = true;
38 
53  // 1-D interpolation is performed
54  // along the x-direction for the u-velocities
55  // along the y-direction for the v-velocities
56 
58  for (int j=0; j < ny; j++)
59  {
60  if(j == 0)
61  {
62  dy0 = 0.5*dy[0];
63  dy1 = 0.5*(dy[0]+dy[1]);
64  }
65  else if(j == ny-1)
66  {
67  dy0 = 0.5*(dy[ny-2]+dy[ny-1]);
68  dy1 = 0.5*dy[ny-1];
69  }
70  else
71  {
72  dy0 = 0.5*(dy[j]+dy[j-1]);
73  dy1 = 0.5*(dy[j]+dy[j+1]);
74  }
75  Cy0 = 2.0 * nu / ( dy1*(dy1+dy0) );
76  Cy1 = 2.0 * nu / ( dy0*(dy1+dy0) );
77 
78  for (int i=0; i < nx-1; i++)
79  {
80  // calculate the row of the matrix
81  I = j*(nx-1) + i;
82 
83  // check if the node is on the boundary or not
84  notBdryNode = (tags[I]==-1);
85 
86  Cx0 = 2.0 * nu / ( dx[i+1]*(dx[i+1]+dx[i]) );
87  Cx1 = 2.0 * nu / ( dx[i]*(dx[i+1]+dx[i]) );
88 
89  // scaling factor (to obtain the normalised matrix)
90  scale = 0.5*(dx[i+1]+dx[i]);
91 
92  // south
93  if(j>0) // no south coefficient for the bottom row
94  {
95  LHost.row_indices[num_elements] = I;
96  LHost.column_indices[num_elements] = I-(nx-1);
97  LHost.values[num_elements] = 0.0;
98 
99  if(notBdryNode) // if the interpolation is definitely not taking place
100  LHost.values[num_elements] = Cy1 * scale / dy[j-1];
101  else if(tags[I]==I-(nx-1) && interpType==LINEAR)
102  LHost.values[num_elements] = coeffs[I] /dt * scale / dy[j-1];
103  else if((tags2[I]==I-2 || tags2[I]==I-2*(nx-1)) && interpType==QUADRATIC)
104  {
105  LHost.column_indices[num_elements] = tags2[I];
106  LHost.values[num_elements] = coeffs2[I] /dt * scale / dy[j-1];
107  }
108 
109  num_elements++;
110  }
111  // west
112  if(i>0) // no west coefficient for the leftmost column
113  {
114  LHost.row_indices[num_elements] = I;
115  LHost.column_indices[num_elements] = I-1;
116  LHost.values[num_elements] = 0.0;
117 
118  if(notBdryNode)
119  LHost.values[num_elements] = Cx1 * scale / dy[j];
120  else if(tags[I]==I-1)
121  LHost.values[num_elements] = coeffs[I] /dt * scale / dy[j];
122  else if(tags[I]==I-(nx-1) && interpType==QUADRATIC)
123  {
124  LHost.column_indices[num_elements] = I-(nx-1);
125  LHost.values[num_elements] = coeffs[I] /dt * scale / dy[j];
126  }
127 
128  num_elements++;
129  }
130  // centre
131  LHost.row_indices[num_elements] = I;
132  LHost.column_indices[num_elements] = I;
133  if(notBdryNode)
134  LHost.values[num_elements] = (-Cy0-Cy1-Cx0-Cx1) * scale / dy[j];
135  else
136  LHost.values[num_elements] = 0.0;
137  num_elements++;
138  // east
139  if(i<nx-2) // no east coefficient for the rightmost column
140  {
141  LHost.row_indices[num_elements] = I;
142  LHost.column_indices[num_elements] = I+1;
143  LHost.values[num_elements] = 0.0;
144 
145  if(notBdryNode)
146  LHost.values[num_elements] = Cx0 * scale / dy[j];
147  else if(tags[I]==I+1)
148  LHost.values[num_elements] = coeffs[I] /dt * scale / dy[j];
149  else if(tags[I]==I+(nx-1) && interpType==QUADRATIC)
150  {
151  LHost.column_indices[num_elements] = I+(nx-1);
152  LHost.values[num_elements] = coeffs[I] /dt * scale / dy[j];
153  }
154 
155  num_elements++;
156  }
157  // north
158  if(j<ny-1)
159  {
160  LHost.row_indices[num_elements] = I;
161  LHost.column_indices[num_elements] = I+(nx-1);
162  LHost.values[num_elements] = 0.0;
163 
164  if(notBdryNode)
165  LHost.values[num_elements] = Cy0 * scale / dy[j+1];
166  else if(tags[I]==I+(nx-1) && interpType==LINEAR)
167  LHost.values[num_elements] = coeffs[I] /dt * scale / dy[j+1];
168  else if((tags2[I]==I+2 || tags2[I]==I+2*(nx-1)) && interpType==QUADRATIC)
169  {
170  LHost.column_indices[num_elements] = tags2[I];
171  LHost.values[num_elements] = coeffs2[I] /dt * scale / dy[j+1];
172  }
173 
174  num_elements++;
175  }
176  }
177  }
178 
180  for (int j=0; j < ny-1; j++)
181  {
182  Cy0 = 2.0 * nu / ( dy[j+1]*(dy[j+1]+dy[j]) );
183  Cy1 = 2.0 * nu / ( dy[j]*(dy[j+1]+dy[j]) );
184 
185  for (int i=0; i < nx; i++)
186  {
187  I = j*nx + i + numU; // calculate the row of the matrix
188  notBdryNode = (tags[I]==-1);
189 
190  if((i>0) && (i<nx-1))
191  {
192  dx0 = 0.5*(dx[i]+dx[i-1]);
193  dx1 = 0.5*(dx[i]+dx[i+1]);
194  }
195  else if(i == 0)
196  {
197  dx0 = 0.5*dx[i];
198  dx1 = 0.5*(dx[i]+dx[i+1]);
199  }
200  else if(i == nx-1)
201  {
202  dx0 = 0.5*(dx[i]+dx[i-1]);
203  dx1 = 0.5*dx[i];
204  }
205  Cx0 = 2.0 * nu /( dx1*(dx1+dx0) );
206  Cx1 = 2.0 * nu /( dx0*(dx1+dx0) );
207 
208  scale = 0.5*(dy[j+1]+dy[j]); // scaling factor (to obtain the normalised matrix)
209 
210  // south
211  if(j>0) // no south coefficient for the bottom row
212  {
213  LHost.row_indices[num_elements] = I;
214  LHost.column_indices[num_elements] = I-nx;
215  LHost.values[num_elements] = 0.0;
216 
217  if(notBdryNode)
218  LHost.values[num_elements] = Cy1 * scale / dx[i];
219  else if(tags[I]==I-nx && interpType==LINEAR)
220  LHost.values[num_elements] = coeffs[I] /dt * scale / dx[i];
221  else if((tags2[I]==I-2 || tags2[I]==I-2*nx) && interpType==QUADRATIC)
222  {
223  LHost.column_indices[num_elements] = tags2[I];
224  LHost.values[num_elements] = coeffs2[I] /dt * scale / dx[i];
225  }
226 
227  num_elements++;
228  }
229  // west
230  if(i>0) // no west coefficient for the leftmost column
231  {
232  LHost.row_indices[num_elements] = I;
233  LHost.column_indices[num_elements] = I-1;
234  LHost.values[num_elements] = 0.0;
235 
236  if(notBdryNode)
237  LHost.values[num_elements] = Cx1 * scale / dx[i-1];
238  else if(tags[I]==I-1)
239  LHost.values[num_elements] = coeffs[I] /dt * scale / dx[i-1];
240  else if(tags[I]==I-nx && interpType==QUADRATIC)
241  {
242  LHost.column_indices[num_elements] = I-nx;
243  LHost.values[num_elements] = coeffs[I] /dt * scale / dx[i-1];
244  }
245 
246  num_elements++;
247  }
248  // centre
249  LHost.row_indices[num_elements] = I;
250  LHost.column_indices[num_elements] = I;
251  if(notBdryNode)
252  LHost.values[num_elements] = (-Cy0-Cy1-Cx0-Cx1) * scale / dx[i];
253  else
254  LHost.values[num_elements] = 0.0;
255  num_elements++;
256  // east
257  if(i<nx-1) // no east coefficient for the rightmost column
258  {
259  LHost.row_indices[num_elements] = I;
260  LHost.column_indices[num_elements] = I+1;
261  LHost.values[num_elements] = 0.0;
262 
263  if(notBdryNode)
264  LHost.values[num_elements] = Cx0 * scale / dx[i+1];
265  else if(tags[I]==I+1)
266  LHost.values[num_elements] = coeffs[I] /dt * scale / dx[i+1];
267  else if(tags[I]==I+nx && interpType==QUADRATIC)
268  {
269  LHost.column_indices[num_elements] = I+nx;
270  LHost.values[num_elements] = coeffs[I] /dt * scale / dx[i+1];
271  }
272 
273  num_elements++;
274  }
275  // north
276  if(j<ny-2)
277  {
278  LHost.row_indices[num_elements] = I;
279  LHost.column_indices[num_elements] = I+nx;
280  LHost.values[num_elements] = 0.0;
281 
282  if(notBdryNode)
283  LHost.values[num_elements] = Cy0 * scale / dx[i];
284  else if(tags[I]==I+nx && interpType==LINEAR)
285  LHost.values[num_elements] = coeffs[I] /dt * scale / dx[i];
286  else if((tags2[I]==I+2 || tags2[I]==I+2*nx) && interpType==QUADRATIC)
287  {
288  LHost.column_indices[num_elements] = tags2[I];
289  LHost.values[num_elements] = coeffs2[I] /dt * scale / dx[i];
290  }
291 
292  num_elements++;
293  }
294  }
295  }
296  L = LHost;
297 #endif
298 } // generateL
299 
300 
301 /*
302 template <typename Matrix, typename Vector>
303 template <>
304 void NavierStokesSolver<Matrix, Vector>::generateBN<3>()
305 {
306  Matrix temp1, temp2;
307  cusp::multiply(Minv, L, temp1);
308  cusp::multiply(temp1, Minv, BN);
309  cusp::add(Minv, BN, BN);
310  cusp::multiply(temp1, BN, temp2);
311  cusp::add(Minv, temp2, BN);
312 }*/
real
double real
Is a float or a double depending on the machine precision.
Definition: types.h:116
QUADRATIC
quadratic
Definition: types.h:96
interpolationType
interpolationType
Specifies the type of interpolation.
Definition: types.h:92
DirectForcingSolver::generateL
virtual void generateL()
LINEAR
linear
Definition: types.h:95
cooH
cusp::coo_matrix< int, real, host_memory > cooH
COO matrix stored on the host.
Definition: types.h:126
Generated on Wed Jul 19 2017 12:27:42 for cuIBM by   doxygen 1.8.11