aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--boltz.cc112
1 files changed, 64 insertions, 48 deletions
diff --git a/boltz.cc b/boltz.cc
index 470d5d8..bf98271 100644
--- a/boltz.cc
+++ b/boltz.cc
@@ -2,6 +2,7 @@
#include <numeric>
#include <iostream>
#include <fstream>
+#include <memory>
#include "src/vector.h"
@@ -48,33 +49,48 @@ struct Cell : DataCell {
}
};
+class CellBuffer {
+ private:
+ const std::size_t dim_x_;
+ const std::size_t dim_y_;
+
+ std::unique_ptr<Cell[]> curr_;
+ std::unique_ptr<Cell[]> prev_;
+
+ public:
+ CellBuffer(std::size_t dimX, std::size_t dimY):
+ dim_x_(dimX),
+ dim_y_(dimY),
+ curr_(new Cell[dimX*dimY]),
+ prev_(new Cell[dimX*dimY]) { }
+
+ void swap() {
+ curr_.swap(prev_);
+ }
+
+ inline Cell& curr(std::size_t x, std::size_t y) {
+ return curr_[y*dim_x_ + x];
+ }
+
+ inline Cell& prev(std::size_t x, std::size_t y) {
+ return prev_[y*dim_x_ + x];
+ }
+
+};
+
constexpr std::size_t dimX = 128;
constexpr std::size_t dimY = 128;
-constexpr double tau = 1.0;
-
-Cell cellsA[dimX][dimY];
-Cell cellsB[dimX][dimY];
+constexpr double tau = 0.6;
-Cell (*newCells)[dimX][dimY] = &cellsA;
-Cell (*oldCells)[dimX][dimY] = &cellsB;
+CellBuffer pop(dimX, dimY);
Density density [dimX][dimY];
Velocity velocity[dimX][dimY];
Force force [dimX][dimY];
-void swapCellBuffers() {
- if ( newCells == &cellsA && oldCells == &cellsB ) {
- newCells = &cellsB;
- oldCells = &cellsA;
- } else {
- newCells = &cellsA;
- oldCells = &cellsB;
- }
-}
-
void computeLbmStep(std::size_t t) {
- swapCellBuffers();
+ pop.swap();
for ( std::size_t x = 1; x < dimX - 1; ++x ) {
for ( std::size_t y = 1; y < dimY - 1; ++y ) {
@@ -85,10 +101,10 @@ void computeLbmStep(std::size_t t) {
Cell eq;
eq.equilibrize(density[x][y], velocity[x][y]);
- // collide & stream
+ // collide (BGK, relax towards equilibrium) & stream
for ( int i = -1; i <= 1; ++i ) {
for ( int j = -1; j <= 1; ++j ) {
- (*newCells)[x+i][y+j].get(i,j) = (*oldCells)[x][y].get(i,j) + 1./tau * (eq.get(i,j) - (*oldCells)[x][y].get(i,j));
+ pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j) + 1./tau * (eq.get(i,j) - pop.prev(x,y).get(i,j));
}
}
}
@@ -96,46 +112,46 @@ void computeLbmStep(std::size_t t) {
// straight wall cell bounce back
for ( std::size_t x = 2; x < dimX - 2; ++x ) {
- (*newCells)[x][1].get( 0, 1) = (*newCells)[x][0].get( 0,-1);
- (*newCells)[x][1].get( 1, 1) = (*newCells)[x][0].get(-1,-1);
- (*newCells)[x][1].get(-1, 1) = (*newCells)[x][0].get( 1,-1);
+ pop.curr(x,1).get( 0, 1) = pop.curr(x,0).get( 0,-1);
+ pop.curr(x,1).get( 1, 1) = pop.curr(x,0).get(-1,-1);
+ pop.curr(x,1).get(-1, 1) = pop.curr(x,0).get( 1,-1);
- (*newCells)[x][dimY-2].get( 0,-1) = (*newCells)[x][dimY-2].get( 0, 1);
- (*newCells)[x][dimY-2].get(-1,-1) = (*newCells)[x][dimY-2].get( 1, 1);
- (*newCells)[x][dimY-2].get( 1,-1) = (*newCells)[x][dimY-2].get(-1, 1);
+ pop.curr(x,dimY-2).get( 0,-1) = pop.curr(x,dimY-2).get( 0, 1);
+ pop.curr(x,dimY-2).get(-1,-1) = pop.curr(x,dimY-2).get( 1, 1);
+ pop.curr(x,dimY-2).get( 1,-1) = pop.curr(x,dimY-2).get(-1, 1);
}
for ( std::size_t y = 2; y < dimY - 2; ++y ) {
- (*newCells)[1][y].get( 1, 0) = (*newCells)[0][y].get(-1, 0);
- (*newCells)[1][y].get( 1,-1) = (*newCells)[0][y].get(-1, 1);
- (*newCells)[1][y].get( 1, 1) = (*newCells)[0][y].get(-1,-1);
+ pop.curr(1,y).get( 1, 0) = pop.curr(0,y).get(-1, 0);
+ pop.curr(1,y).get( 1,-1) = pop.curr(0,y).get(-1, 1);
+ pop.curr(1,y).get( 1, 1) = pop.curr(0,y).get(-1,-1);
- (*newCells)[dimX-2][y].get(-1, 0) = (*newCells)[dimX - 1][y].get( 1, 0);
- (*newCells)[dimX-2][y].get(-1, 1) = (*newCells)[dimX - 1][y].get( 1,-1);
- (*newCells)[dimX-2][y].get(-1,-1) = (*newCells)[dimX - 1][y].get( 1, 1);
+ pop.curr(dimX-2,y).get(-1, 0) = pop.curr(dimX-1,y).get( 1, 0);
+ pop.curr(dimX-2,y).get(-1, 1) = pop.curr(dimX-1,y).get( 1,-1);
+ pop.curr(dimX-2,y).get(-1,-1) = pop.curr(dimX-1,y).get( 1, 1);
}
// edge wall cell bounce back
- (*newCells)[1][1].get( 0, 1) = (*newCells)[0][0].get( 0,-1);
- (*newCells)[1][1].get( 1, 1) = (*newCells)[0][0].get(-1,-1);
- (*newCells)[1][1].get( 1, 0) = (*newCells)[0][0].get(-1, 0);
+ pop.curr(1,1).get( 0, 1) = pop.curr(0,0).get( 0,-1);
+ pop.curr(1,1).get( 1, 1) = pop.curr(0,0).get(-1,-1);
+ pop.curr(1,1).get( 1, 0) = pop.curr(0,0).get(-1, 0);
- (*newCells)[1][dimY-2].get( 0,-1) = (*newCells)[0][dimY-1].get( 0, 1);
- (*newCells)[1][dimY-2].get( 1,-1) = (*newCells)[0][dimY-1].get(-1, 1);
- (*newCells)[1][dimY-2].get( 1, 0) = (*newCells)[0][dimY-1].get(-1, 0);
+ pop.curr(1,dimY-2).get( 0,-1) = pop.curr(0,dimY-1).get( 0, 1);
+ pop.curr(1,dimY-2).get( 1,-1) = pop.curr(0,dimY-1).get(-1, 1);
+ pop.curr(1,dimY-2).get( 1, 0) = pop.curr(0,dimY-1).get(-1, 0);
- (*newCells)[dimX-2][1].get( 0, 1) = (*newCells)[dimX-1][0].get( 0,-1);
- (*newCells)[dimX-2][1].get(-1, 1) = (*newCells)[dimX-1][0].get( 1,-1);
- (*newCells)[dimX-2][1].get(-1, 0) = (*newCells)[dimX-1][0].get( 1, 0);
+ pop.curr(dimX-2,1).get( 0, 1) = pop.curr(dimX-1,0).get( 0,-1);
+ pop.curr(dimX-2,1).get(-1, 1) = pop.curr(dimX-1,0).get( 1,-1);
+ pop.curr(dimX-2,1).get(-1, 0) = pop.curr(dimX-1,0).get( 1, 0);
- (*newCells)[dimX-2][dimY-2].get( 0,-1) = (*newCells)[dimX-1][dimY-1].get( 0, 1);
- (*newCells)[dimX-2][dimY-2].get(-1,-1) = (*newCells)[dimX-1][dimY-1].get( 1, 1);
- (*newCells)[dimX-2][dimY-2].get(-1, 0) = (*newCells)[dimX-1][dimY-1].get( 1, 0);
+ pop.curr(dimX-2,dimY-2).get( 0,-1) = pop.curr(dimX-1,dimY-1).get( 0, 1);
+ pop.curr(dimX-2,dimY-2).get(-1,-1) = pop.curr(dimX-1,dimY-1).get( 1, 1);
+ pop.curr(dimX-2,dimY-2).get(-1, 0) = pop.curr(dimX-1,dimY-1).get( 1, 0);
// update density, velocity field
for ( std::size_t x = 0; x < dimX; ++x ) {
for ( std::size_t y = 0; y < dimY; ++y ) {
- density[x][y] = (*newCells)[x][y].sum();
- velocity[x][y] = (*newCells)[x][y].velocity(density[x][y]);
+ density[x][y] = pop.curr(x,y).sum();
+ velocity[x][y] = pop.curr(x,y).velocity(density[x][y]);
}
}
}
@@ -189,8 +205,8 @@ void init() {
velocity[x][y] = { 0.0, 0.0 };
force[x][y] = { 0.0, 0.0 };
- (*newCells)[x][y].equilibrize(density[x][y], velocity[x][y]);
- (*oldCells)[x][y].equilibrize(density[x][y], velocity[x][y]);
+ pop.curr(x,y).equilibrize(density[x][y], velocity[x][y]);
+ pop.prev(x,y).equilibrize(density[x][y], velocity[x][y]);
}
}
@@ -207,7 +223,7 @@ void init() {
int main() {
init();
- for ( std::size_t t = 0; t < 100; ++t ) {
+ for ( std::size_t t = 0; t < 500; ++t ) {
computeLbmStep(t);
if ( t % 2 == 0 ) {