cuIBM
A GPU-based Immersed Boundary Method code
DFImprovedSolver.cu
Go to the documentation of this file.
1 
7 #include <cusp/io/matrix_market.h>
8 #include <cusp/blas/blas.h>
9 
10 #include "DFImprovedSolver.h"
11 #include "solvers/kernels/generateQT.h"
12 
13 
17 template <typename memoryType>
18 DFImprovedSolver<memoryType>::DFImprovedSolver(parameterDB *pDB, domain *dInfo)
19 {
20  NavierStokesSolver<memoryType>::paramDB = pDB;
21  NavierStokesSolver<memoryType>::domInfo = dInfo;
22 } // DFImprovedSolver
23 
24 
28 template <typename memoryType>
29 void DFImprovedSolver<memoryType>::generateQT()
30 {
31  int nx = NavierStokesSolver<memoryType>::domInfo->nx,
32  ny = NavierStokesSolver<memoryType>::domInfo->ny;
33 
34  const int N_u = (nx-1)*ny;
35 
36  NavierStokesSolver<memoryType>::generateQT();
37  DirectForcingSolver<memoryType>::updateQ();
38 
39  cusp::coo_matrix<int, real, host_memory> QTHost(nx*ny, (nx-1)*ny+nx*(ny-1), 4*nx*ny-2*(nx+ny));
40  cusp::blas::fill(QTHost.row_indices, -1);
41  cusp::blas::fill(QTHost.column_indices, -1);
42  cusp::blas::fill(QTHost.values, 0.0);
43 
44  int idx = 0;
45  for(int j=0; j<ny; j++)
46  {
47  for(int i=0; i<nx; i++)
48  {
49  int row = j*nx+i;
50  if(i>0)
51  {
52  int I = j*(nx-1)+(i-1);
53  if(DirectForcingSolver<memoryType>::tags[I]==-1)
54  {
55  QTHost.row_indices[idx] = row;
56  QTHost.column_indices[idx] = I;
57  QTHost.values[idx] = 1.0;
58  idx++;
59  }
60  else
61  {
62  bool flag = false;
63  int start;
64  start = (idx>4)? idx-4 : 0;
65  for(int l=start; l<idx && !flag; l++)
66  {
67  if(QTHost.row_indices[l]==row && QTHost.column_indices[l]==DirectForcingSolver<memoryType>::tags[I])
68  {
69  flag = true;
70  QTHost.values[l] += DirectForcingSolver<memoryType>::coeffs[I];
71  }
72  }
73  if(!flag)
74  {
75  QTHost.row_indices[idx] = row;
76  QTHost.column_indices[idx] = DirectForcingSolver<memoryType>::tags[I];
77  QTHost.values[idx] = DirectForcingSolver<memoryType>::coeffs[I];
78  idx++;
79  }
80  }
81  }
82  if(i<nx-1)
83  {
84  int I = j*(nx-1)+i;
85  if(DirectForcingSolver<memoryType>::tags[I]==-1)
86  {
87  QTHost.row_indices[idx] = row;
88  QTHost.column_indices[idx] = I;
89  QTHost.values[idx] = -1.0;
90  idx++;
91  }
92  else
93  {
94  bool flag = false;
95  int start;
96  start = (idx>4)? idx-4 : 0;
97  for(int l=start; l<idx && !flag; l++)
98  {
99  if(QTHost.row_indices[l]==row && QTHost.column_indices[l]==DirectForcingSolver<memoryType>::tags[I])
100  {
101  flag = true;
102  QTHost.values[l] -= DirectForcingSolver<memoryType>::coeffs[I];
103  }
104  }
105  if(!flag)
106  {
107  QTHost.row_indices[idx] = row;
108  QTHost.column_indices[idx] = DirectForcingSolver<memoryType>::tags[I];
109  QTHost.values[idx] = -DirectForcingSolver<memoryType>::coeffs[I];
110  idx++;
111  }
112  }
113  }
114  if(j>0)
115  {
116  int I = (j-1)*nx+i+N_u;
117  if(DirectForcingSolver<memoryType>::tags[I]==-1)
118  {
119  QTHost.row_indices[idx] = row;
120  QTHost.column_indices[idx] = I;
121  QTHost.values[idx] = 1.0;
122  idx++;
123  }
124  else
125  {
126  bool flag = false;
127  int start;
128  start = (idx>4)? idx-4 : 0;
129  for(int l=start; l<idx && !flag; l++)
130  {
131  if(QTHost.row_indices[l]==row && QTHost.column_indices[l]==DirectForcingSolver<memoryType>::tags[I])
132  {
133  flag = true;
134  QTHost.values[l] += DirectForcingSolver<memoryType>::coeffs[I];
135  }
136  }
137  if(!flag)
138  {
139  QTHost.row_indices[idx] = row;
140  QTHost.column_indices[idx] = DirectForcingSolver<memoryType>::tags[I];
141  QTHost.values[idx] = DirectForcingSolver<memoryType>::coeffs[I];
142  idx++;
143  }
144  }
145  }
146  if(j<ny-1)
147  {
148  int I = j*nx+i+N_u;
149  if(DirectForcingSolver<memoryType>::tags[I]==-1)
150  {
151  QTHost.row_indices[idx] = row;
152  QTHost.column_indices[idx] = I;
153  QTHost.values[idx] = -1.0;
154  idx++;
155  }
156  else
157  {
158  bool flag = false;
159  int start;
160  start = (idx>4)? idx-4 : 0;
161  for(int l=start; l<idx && !flag; l++)
162  {
163  if(QTHost.row_indices[l]==row && QTHost.column_indices[l]==DirectForcingSolver<memoryType>::tags[I])
164  {
165  flag = true;
166  QTHost.values[l] -= DirectForcingSolver<memoryType>::coeffs[I];
167  }
168  }
169  if(!flag)
170  {
171  QTHost.row_indices[idx] = row;
172  QTHost.column_indices[idx] = DirectForcingSolver<memoryType>::tags[I];
173  QTHost.values[idx] = -DirectForcingSolver<memoryType>::coeffs[I];
174  idx++;
175  }
176  }
177  }
178  }
179  }
180  NavierStokesSolver<memoryType>::QT = QTHost;
181  NavierStokesSolver<memoryType>::QT.resize(nx*ny, (nx-1)*ny+nx*(ny-1), idx);
182  /*std::cout << "\nQT stuff:\n";
183  std::cout << "Copied and resized matrix." << std::endl;
184  std::cout << "Original size: " << QTHost.values.size() << std::endl;
185  std::cout << "Actual size : " << idx << std::endl;
186  cusp::io::write_matrix_market_file(NavierStokesSolver<memoryType>::Q, "Q.mtx");
187  std::cout << "Wrote Q to file." << std::endl;
188  cusp::io::write_matrix_market_file(NavierStokesSolver<memoryType>::QT, "QT.mtx");
189  std::cout << "Wrote QT to file." << std::endl;*/
190 } // generateQT
191 
192 
193 /*
194 template <typename memoryType>
195 void DFImprovedSolver<memoryType>::generateC()
196 {
197  int nx = NavierStokesSolver<memoryType>::domInfo->nx,
198  ny = NavierStokesSolver<memoryType>::domInfo->ny;
199 
200  parameterDB &db = *NavierStokesSolver<memoryType>::paramDB;
201  real dt = db["simulation"]["dt"].get<real>();
202 
203  const int ii = 2, jj = 2;
204  const int N_u = (nx-1)*ny;
205 
206  int isColumnNonZero[5][5];
207  for(int m=-2; m<=2; m++)
208  {
209  for(int l=-2; l<=2; l++)
210  {
211  isColumnNonZero[jj+m][ii+l]=0;
212  }
213  }
214 
215  int num_nonzeros = 0;
216  for(int j=0; j<ny; j++)
217  {
218  for(int i=0; i<nx; i++)
219  {
220  if(j>0)
221  {
222  int I = (j-1)*nx+i+N_u;
223  if(DirectForcingSolver<memoryType>::tags[I]==-1)
224  {
225  isColumnNonZero[jj-1][ii] += 1;
226  isColumnNonZero[jj][ii] += 1;
227  }
228  else
229  {
230  int diff = DirectForcingSolver<memoryType>::tags[I]-I;
231  if(diff < -1)
232  {
233  isColumnNonZero[jj-2][ii] += 1;
234  isColumnNonZero[jj-1][ii] += 1;
235  }
236  else if(diff == -1)
237  {
238  isColumnNonZero[jj-1][ii-1] += 1;
239  isColumnNonZero[jj][ii-1] += 1;
240  }
241  else if(diff == 1)
242  {
243  isColumnNonZero[jj-1][ii+1] += 1;
244  isColumnNonZero[jj][ii+1] += 1;
245  }
246  else if(diff > 1)
247  {
248  isColumnNonZero[jj][ii] += 1;
249  isColumnNonZero[jj+1][ii] += 1;
250  }
251  }
252  }
253  if(i>0)
254  {
255  int I = j*(nx-1)+(i-1);
256  if(DirectForcingSolver<memoryType>::tags[I]==-1)
257  {
258  isColumnNonZero[jj][ii-1] += 1;
259  isColumnNonZero[jj][ii] += 1;
260  }
261  else
262  {
263  int diff = DirectForcingSolver<memoryType>::tags[I]-I;
264  if(diff < -1)
265  {
266  isColumnNonZero[jj-1][ii-1] += 1;
267  isColumnNonZero[jj-1][ii] += 1;
268  }
269  else if(diff == -1)
270  {
271  isColumnNonZero[jj][ii-2] += 1;
272  isColumnNonZero[jj][ii-1] += 1;
273  }
274  else if(diff == 1)
275  {
276  isColumnNonZero[jj][ii] += 1;
277  isColumnNonZero[jj][ii+1] += 1;
278  }
279  else if(diff > 1)
280  {
281  isColumnNonZero[jj+1][ii-1] += 1;
282  isColumnNonZero[jj+1][ii] += 1;
283  }
284  }
285  }
286  if(i<nx-1)
287  {
288  int I = j*(nx-1)+i;
289  if(DirectForcingSolver<memoryType>::tags[I]==-1)
290  {
291  isColumnNonZero[jj][ii+1] += 1;
292  isColumnNonZero[jj][ii] += 1;
293  }
294  else
295  {
296  int diff = DirectForcingSolver<memoryType>::tags[I]-I;
297  if(diff < -1)
298  {
299  isColumnNonZero[jj-1][ii+1] += 1;
300  isColumnNonZero[jj-1][ii] += 1;
301  }
302  else if(diff == -1)
303  {
304  isColumnNonZero[jj][ii] += 1;
305  isColumnNonZero[jj][ii-1] += 1;
306  }
307  else if(diff == 1)
308  {
309  isColumnNonZero[jj][ii+2] += 1;
310  isColumnNonZero[jj][ii+1] += 1;
311  }
312  else if(diff > 1)
313  {
314  isColumnNonZero[jj+1][ii+1] += 1;
315  isColumnNonZero[jj+1][ii] += 1;
316  }
317  }
318  }
319  if(j<ny-1)
320  {
321  int I = j*nx+i+N_u;
322  if(DirectForcingSolver<memoryType>::tags[I]==-1)
323  {
324  isColumnNonZero[jj+1][ii] += 1;
325  isColumnNonZero[jj][ii] += 1;
326  }
327  else
328  {
329  int diff = DirectForcingSolver<memoryType>::tags[I]-I;
330  if(diff < -1)
331  {
332  isColumnNonZero[jj][ii] += 1;
333  isColumnNonZero[jj-1][ii] += 1;
334  }
335  else if(diff == -1)
336  {
337  isColumnNonZero[jj+1][ii-1] += 1;
338  isColumnNonZero[jj][ii-1] += 1;
339  }
340  else if(diff == 1)
341  {
342  isColumnNonZero[jj+1][ii+1] += 1;
343  isColumnNonZero[jj][ii+1] += 1;
344  }
345  else if(diff > 1)
346  {
347  isColumnNonZero[jj+2][ii] += 1;
348  isColumnNonZero[jj+1][ii] += 1;
349  }
350  }
351  }
352  int numNonZeroColumns = 0;
353  //std::cout << "(" << i << "," << j << ")|";
354  for(int m=-2; m<=2; m++)
355  {
356  for(int l=-2; l<=2; l++)
357  {
358  //std::cout << isColumnNonZero[jj+m][ii+l] << ",";
359  if(isColumnNonZero[jj+m][ii+l]) numNonZeroColumns++;
360  isColumnNonZero[jj+m][ii+l] = 0;
361  }
362  //std::cout << "|";
363  }
364  //std::cout << numNonZeroColumns << std::endl;
365  num_nonzeros += numNonZeroColumns;
366  }
367  }
368  //std::cout << "Total nonzeros: " << num_nonzeros << std::endl;
369 
370  cusp::coo_matrix<int, real, host_memory> CHost(nx*ny, nx*ny, num_nonzeros);
371 
372  real valuesInColumns[5][5];
373 
374  for(int m=-2; m<=2; m++)
375  {
376  for(int l=-2; l<=2; l++)
377  {
378  valuesInColumns[jj+m][ii+l]=0.0;
379  }
380  }
381 
382  int idx = 0;
383  for(int j=0; j<ny; j++)
384  {
385  for(int i=0; i<nx; i++)
386  {
387  if(j>0)
388  {
389  int I = (j-1)*nx+i+N_u;
390  if(DirectForcingSolver<memoryType>::tags[I]==-1)
391  {
392  isColumnNonZero[jj-1][ii] += 1;
393  valuesInColumns[jj-1][ii] -= 1.0;
394  isColumnNonZero[jj][ii] += 1;
395  valuesInColumns[jj][ii] += 1.0;
396  }
397  else
398  {
399  int diff = DirectForcingSolver<memoryType>::tags[I]-I;
400  real xi = DirectForcingSolver<memoryType>::coeffs[I];
401  if(diff < -1)
402  {
403  isColumnNonZero[jj-2][ii] += 1;
404  valuesInColumns[jj-2][ii] -= xi;
405  isColumnNonZero[jj-1][ii] += 1;
406  valuesInColumns[jj-1][ii] += xi;
407  }
408  else if(diff == -1)
409  {
410  isColumnNonZero[jj-1][ii-1] += 1;
411  valuesInColumns[jj-1][ii-1] -= xi;
412  isColumnNonZero[jj][ii-1] += 1;
413  valuesInColumns[jj][ii-1] += xi;
414  }
415  else if(diff == 1)
416  {
417  isColumnNonZero[jj-1][ii+1] += 1;
418  valuesInColumns[jj-1][ii+1] -= xi;
419  isColumnNonZero[jj][ii+1] += 1;
420  valuesInColumns[jj][ii+1] += xi;
421  }
422  else if(diff > 1)
423  {
424  isColumnNonZero[jj][ii] += 1;
425  valuesInColumns[jj][ii] -= xi;
426  isColumnNonZero[jj+1][ii] += 1;
427  valuesInColumns[jj+1][ii] += xi;
428  }
429  }
430  }
431  if(i>0)
432  {
433  int I = j*(nx-1)+(i-1);
434  if(DirectForcingSolver<memoryType>::tags[I]==-1)
435  {
436  isColumnNonZero[jj][ii-1] += 1;
437  valuesInColumns[jj][ii-1] -= 1.0;
438  isColumnNonZero[jj][ii] += 1;
439  valuesInColumns[jj][ii] += 1.0;
440  }
441  else
442  {
443  int diff = DirectForcingSolver<memoryType>::tags[I]-I;
444  real xi = DirectForcingSolver<memoryType>::coeffs[I];
445  if(diff < -1)
446  {
447  isColumnNonZero[jj-1][ii-1] += 1;
448  valuesInColumns[jj-1][ii-1] -= xi;
449  isColumnNonZero[jj-1][ii] += 1;
450  valuesInColumns[jj-1][ii] += xi;
451  }
452  else if(diff == -1)
453  {
454  isColumnNonZero[jj][ii-2] += 1;
455  valuesInColumns[jj][ii-2] -= xi;
456  isColumnNonZero[jj][ii-1] += 1;
457  valuesInColumns[jj][ii-1] += xi;
458  }
459  else if(diff == 1)
460  {
461  isColumnNonZero[jj][ii] += 1;
462  valuesInColumns[jj][ii] -= xi;
463  isColumnNonZero[jj][ii+1] += 1;
464  valuesInColumns[jj][ii+1] += xi;
465  }
466  else if(diff > 1)
467  {
468  isColumnNonZero[jj+1][ii-1] += 1;
469  valuesInColumns[jj+1][ii-1] -= xi;
470  isColumnNonZero[jj+1][ii] += 1;
471  valuesInColumns[jj+1][ii] += xi;
472  }
473  }
474  }
475  if(i<nx-1)
476  {
477  int I = j*(nx-1)+i;
478  if(DirectForcingSolver<memoryType>::tags[I]==-1)
479  {
480  isColumnNonZero[jj][ii+1] += 1;
481  valuesInColumns[jj][ii+1] -= 1.0;
482  isColumnNonZero[jj][ii] += 1;
483  valuesInColumns[jj][ii] += 1.0;
484  }
485  else
486  {
487  int diff = DirectForcingSolver<memoryType>::tags[I]-I;
488  real xi = DirectForcingSolver<memoryType>::coeffs[I];
489  if(diff < -1)
490  {
491  isColumnNonZero[jj-1][ii+1] += 1;
492  valuesInColumns[jj-1][ii+1] -= xi;
493  isColumnNonZero[jj-1][ii] += 1;
494  valuesInColumns[jj-1][ii] += xi;
495  }
496  else if(diff == -1)
497  {
498  isColumnNonZero[jj][ii] += 1;
499  valuesInColumns[jj][ii] -= xi;
500  isColumnNonZero[jj][ii-1] += 1;
501  valuesInColumns[jj][ii-1] += xi;
502  }
503  else if(diff == 1)
504  {
505  isColumnNonZero[jj][ii+2] += 1;
506  valuesInColumns[jj][ii+2] -= xi;
507  isColumnNonZero[jj][ii+1] += 1;
508  valuesInColumns[jj][ii+1] += xi;
509  }
510  else if(diff > 1)
511  {
512  isColumnNonZero[jj+1][ii+1] += 1;
513  valuesInColumns[jj+1][ii+1] -= xi;
514  isColumnNonZero[jj+1][ii] += 1;
515  valuesInColumns[jj+1][ii] += xi;
516  }
517  }
518  }
519  if(j<ny-1)
520  {
521  int I = j*nx+i+N_u;
522  if(DirectForcingSolver<memoryType>::tags[I]==-1)
523  {
524  isColumnNonZero[jj+1][ii] += 1;
525  valuesInColumns[jj+1][ii] -= 1.0;
526  isColumnNonZero[jj][ii] += 1;
527  valuesInColumns[jj][ii] += 1.0;
528  }
529  else
530  {
531  int diff = DirectForcingSolver<memoryType>::tags[I]-I;
532  real xi = DirectForcingSolver<memoryType>::coeffs[I];
533  if(diff < -1)
534  {
535  isColumnNonZero[jj][ii] += 1;
536  valuesInColumns[jj][ii] -= xi;
537  isColumnNonZero[jj-1][ii] += 1;
538  valuesInColumns[jj-1][ii] += xi;
539  }
540  else if(diff == -1)
541  {
542  isColumnNonZero[jj+1][ii-1] += 1;
543  valuesInColumns[jj+1][ii-1] -= xi;
544  isColumnNonZero[jj][ii-1] += 1;
545  valuesInColumns[jj][ii-1] += xi;
546  }
547  else if(diff == 1)
548  {
549  isColumnNonZero[jj+1][ii+1] += 1;
550  valuesInColumns[jj+1][ii+1] -= xi;
551  isColumnNonZero[jj][ii+1] += 1;
552  valuesInColumns[jj][ii+1] += xi;
553  }
554  else if(diff > 1)
555  {
556  isColumnNonZero[jj+2][ii] += 1;
557  valuesInColumns[jj+2][ii] -= xi;
558  isColumnNonZero[jj+1][ii] += 1;
559  valuesInColumns[jj+1][ii] += xi;
560  }
561  }
562  }
563  int row = j*nx+i;
564  for(int m=-2; m<=2; m++)
565  {
566  for(int l=-2; l<=2; l++)
567  {
568  if(isColumnNonZero[jj+m][ii+l])
569  {
570  CHost.row_indices[idx] = row;
571  CHost.column_indices[idx] = row + m*nx + l;
572  CHost.values[idx] = valuesInColumns[jj+m][ii+l];
573  if(CHost.row_indices[idx]==(ny/2)*nx+nx/2 && CHost.row_indices[idx]==CHost.column_indices[idx])
574  {
575  CHost.values[idx]+=CHost.values[idx];
576  }
577  idx++;
578  }
579  isColumnNonZero[jj+m][ii+l] = 0;
580  valuesInColumns[jj+m][ii+l] = 0.0;
581  }
582  }
583  }
584  }
585  CHost.sort_by_row_and_column();
586  CHost.values[0] += CHost.values[0];
587  NavierStokesSolver<memoryType>::C = CHost;
588  //cusp::io::write_matrix_market_file(NavierStokesSolver<memoryType>::C, "C-generateC.mtx");
589  cusp::blas::scal(NavierStokesSolver<memoryType>::C.values, dt);
590 } // generateC
591 */
592 
593 // specialization of the class
594 template class DFImprovedSolver<device_memory>;
DFImprovedSolver
Second-order fully-discrete direct forcing method.
Definition: DFImprovedSolver.h:17
DirectForcingSolver
A fully discrete formulation of the direct forcing method.
Definition: DirectForcingSolver.h:27
parameterDB
std::map< std::string, componentParameter > parameterDB
Map from a string to a componentParameter.
Definition: parameterDB.h:64
domain
Stores information about the computational grid.
Definition: domain.h:16
DFImprovedSolver::generateQT
virtual void generateQT()
Compute the modified divergence operator.
Definition: DFImprovedSolver.cu:29
DFImprovedSolver.h
Definition of the class DFImprovedSolver.
DFImprovedSolver::DFImprovedSolver
DFImprovedSolver(parameterDB *pDB=NULL, domain *dInfo=NULL)
Constructor – get simulation parameters and grid.
Definition: DFImprovedSolver.cu:18
NavierStokesSolver
Solves the Navier-Stokes equations in a rectangular domain.
Definition: NavierStokesSolver.h:26
generateQT.h
Declaration of the kernels to generate gradient matrix and interpolation matrix.
NavierStokesSolver::generateQT
virtual void generateQT()
Generated on Wed Jul 19 2017 12:27:42 for cuIBM by   doxygen 1.8.11