1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
|
#pragma once
#include "data_cell.h"
#include "data_cell_buffer.h"
#include "fluid_buffer.h"
constexpr DataCell weight{
1./36., 1./9., 1./36.,
1./9., 4./9., 1./9.,
1./36., 1./9., 1./36
};
struct Cell : DataCell {
void equilibrize(Density d, Velocity v) {
for ( int i = -1; i <= 1; ++i ) {
for ( int j = -1; j <= 1; ++j ) {
get(i,j) = weight.get(i,j) * d * (1 + 3*v.comp(i,j) + 4.5*sq(v.comp(i,j)) - 1.5*sq(v.norm()));
}
}
}
double sum() const {
return get(-1, 1) + get( 0, 1) + get( 1, 1) + get(-1, 0) + get( 0, 0) + get( 1, 0) + get(-1,-1) + get( 0,-1) + get( 1,-1);
}
Velocity velocity(Density d) const {
return 1./d * Velocity{
get( 1, 0) - get(-1, 0) + get( 1, 1) - get(-1,-1) + get( 1,-1) - get(-1,1),
get( 0, 1) - get( 0,-1) + get( 1, 1) - get(-1,-1) - get( 1,-1) + get(-1,1)
};
}
};
void streamFluidCell(DataCellBuffer& pop, std::size_t x, std::size_t y) {
if ( x != 0 && x != pop.dimX()-1 && y != 0 && y != pop.dimY()-1 ) {
// stream internal cells
for ( int i = -1; i <= 1; ++i ) {
for ( int j = -1; j <= 1; ++j ) {
pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j);
}
}
} else {
// stream boundary cells,
// missing populations to be determined by boundary conditions
for ( int i = -1; i <= 1; ++i ) {
for ( int j = -1; j <= 1; ++j ) {
if ( (x > 0 || i >= 0) && x+i <= pop.dimX()-1 && (y > 0 || j >= 0) && y+j <= pop.dimY()-1 ) {
pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j);
}
}
}
}
}
void collideFluidCell(
double omega,
DataCellBuffer& pop, FluidBuffer& fluid,
std::size_t x, std::size_t y) {
// compute equilibrium
Cell eq;
eq.equilibrize(fluid.density(x,y), fluid.velocity(x,y));
// collide (BGK, relax towards equilibrium)
if ( x != 0 && x != pop.dimX()-1 && y != 0 && y != pop.dimY()-1 ) {
for ( int i = -1; i <= 1; ++i ) {
for ( int j = -1; j <= 1; ++j ) {
pop.curr(x,y).get(i,j) = pop.curr(x,y).get(i,j) + omega * (eq.get(i,j) - pop.curr(x,y).get(i,j));
}
}
} else {
// partial collide for boundary cells
for ( int i = -1; i <= 1; ++i ) {
for ( int j = -1; j <= 1; ++j ) {
if ( (x > 0 || i >= 0) && x+i <= pop.dimX()-1 && (y > 0 || j >= 0) && y+j <= pop.dimY()-1 ) {
pop.curr(x,y).get(i,j) = pop.curr(x,y).get(i,j) + omega * (eq.get(i,j) - pop.curr(x,y).get(i,j));
}
}
}
}
}
|
