aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--boltz.cc133
1 files changed, 61 insertions, 72 deletions
diff --git a/boltz.cc b/boltz.cc
index 29f8518..3796b87 100644
--- a/boltz.cc
+++ b/boltz.cc
@@ -27,7 +27,6 @@ struct DataCell {
};
using Velocity = Vector<double>;
-using Force = Vector<double>;
using Density = double;
constexpr DataCell weight{
@@ -93,13 +92,12 @@ class CellBuffer {
}
};
-constexpr std::size_t dimX = 128;
-constexpr std::size_t dimY = 128;
+constexpr std::size_t dimX = 256;
+constexpr std::size_t dimY = 256;
-constexpr double lidVelocityX = 0.05;
-constexpr double reynolds = 1000;
-constexpr double nu = lidVelocityX * dimX / reynolds;
-constexpr double tau = 3*nu + 0.5;
+constexpr double tau = 2./3.;
+constexpr int reynolds = 1000;
+constexpr double lidVelocityX = reynolds * (2.*tau - 1.) / (6.*(dimX - 1.));
constexpr double omega = 1. / tau;
@@ -107,7 +105,6 @@ CellBuffer pop(dimX, dimY);
Density density [dimX][dimY];
Velocity velocity[dimX][dimY];
-Force force [dimX][dimY];
void computeFluidCell(std::size_t x, std::size_t y) {
// compute equilibrium
@@ -115,9 +112,21 @@ void computeFluidCell(std::size_t x, std::size_t y) {
eq.equilibrize(density[x][y], velocity[x][y]);
// collide (BGK, relax towards equilibrium) & stream
- 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) + omega * (eq.get(i,j) - pop.prev(x,y).get(i,j));
+ if ( x != 0 && x != dimX - 1 && y != 0 && y != dimY -1 ) {
+ 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) + omega * (eq.get(i,j) - pop.prev(x,y).get(i,j));
+ }
+ }
+ } else {
+ // partial collide and stream for boundary cells,
+ // remaining populations to be determined by boundary conditions
+ for ( int i = -1; i <= 1; ++i ) {
+ for ( int j = -1; j <= 1; ++j ) {
+ if ( x+i >= 0 && x+i <= dimX - 1 && y+j >= 0 && y+j <= dimY - 1 ) {
+ pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j) + omega * (eq.get(i,j) - pop.prev(x,y).get(i,j));
+ }
+ }
}
}
}
@@ -144,78 +153,59 @@ std::pair<Vector<int>, Vector<int>> neighbors(Vector<int> v) {
void computeWallCell(Vector<std::size_t> cell, Vector<int> normal) {
const auto [neighborA, neighborB] = neighbors(normal);
- pop.curr(cell).get(neighborA) = pop.curr(cell-neighborA).get(-neighborA);
- pop.curr(cell).get(normal ) = pop.curr(cell-normal ).get(-normal );
- pop.curr(cell).get(neighborB) = pop.curr(cell-neighborB).get(-neighborB);
+ pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA);
+ pop.curr(cell).get(normal ) = pop.curr(cell).get(-normal );
+ pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB);
}
void computeSimpleVelocityWallCell(Vector<std::size_t> cell, Vector<int> normal, double vX) {
const auto [neighborA, neighborB] = neighbors(normal);
- pop.curr(cell).get(neighborA) = pop.curr(cell-neighborA).get(-neighborA) - 6 * weight.get( 1, 1) * density[cell[0]][cell[1]] * vX;
- pop.curr(cell).get(normal ) = pop.curr(cell-normal ).get(-normal );
- pop.curr(cell).get(neighborB) = pop.curr(cell-neighborB).get(-neighborB) + 6 * weight.get(-1, 1) * density[cell[0]][cell[0]] * vX;
+ pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA) - 6 * weight.get( 1, 1) * density[cell[0]][cell[1]] * vX;
+ pop.curr(cell).get(normal ) = pop.curr(cell).get(-normal );
+ pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB) + 6 * weight.get(-1, 1) * density[cell[0]][cell[0]] * vX;
}
void computeZouHeVelocityWallCell(Vector<std::size_t> cell, Vector<int> normal, double vX) {
const auto [neighborA, neighborB] = neighbors(normal);
- const double rho = pop.curr(cell-neighborA).get(-1,0) + pop.curr(cell-normal).get(0,0) + pop.curr(cell-neighborB).get(1,0)
+ const double rho = pop.curr(cell).get(-1,0) + pop.curr(cell).get(0,0) + pop.curr(cell).get(1,0)
+ 2.*(
- pop.curr(cell-neighborA).get(-neighborA) +
- pop.curr(cell-normal ).get(-normal ) +
- pop.curr(cell-neighborB).get(-neighborB)
+ pop.curr(cell).get(-neighborA) +
+ pop.curr(cell).get(-normal ) +
+ pop.curr(cell).get(-neighborB)
);
- const double hor = pop.curr(cell-neighborA).get(-1,0) - pop.curr(cell-neighborB).get(1,0);
- pop.curr(cell).get(neighborA) = pop.curr(cell-neighborA).get(-neighborA) + 0.5*hor - 0.5*vX*rho;
- pop.curr(cell).get(normal ) = pop.curr(cell-normal ).get(-normal );
- pop.curr(cell).get(neighborB) = pop.curr(cell-neighborB).get(-neighborB) - 0.5*hor + 0.5*vX*rho;
+ pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA) + 0.5*(pop.curr(cell).get( 1,0) - pop.curr(cell).get(-1,0) - vX*rho);
+ pop.curr(cell).get(normal ) = pop.curr(cell).get(-normal );
+ pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB) + 0.5*(pop.curr(cell).get(-1,0) - pop.curr(cell).get( 1,0) + vX*rho);
}
void computeLbmStep(std::size_t t) {
pop.swap();
- for ( std::size_t x = 1; x < dimX - 1; ++x ) {
- for ( std::size_t y = 1; y < dimY - 1; ++y ) {
- //if ( x <= 20 || x >= 40 || y <= 20 || y >= 80 ) {
- computeFluidCell(x, y);
- //}
+ for ( std::size_t x = 1; x < dimX-1; ++x ) {
+ for ( std::size_t y = 0; y < dimY; ++y ) {
+ computeFluidCell(x, y);
}
}
- // obstacle
- /*{
- for ( std::size_t x = 21; x < 40; ++x ) {
- computeWallCell(x, 80, { 0, 1});
- computeWallCell(x, 20, { 0,-1});
- }
- for ( std::size_t y = 21; y < 80; ++y ) {
- computeWallCell(40, y, { 1, 0});
- computeWallCell(20, y, {-1, 0});
- }
-
- computeWallCell(40, 80, { 1, 1});
- computeWallCell(40, 20, { 1,-1});
- computeWallCell(20, 80, {-1, 1});
- computeWallCell(20, 20, {-1,-1});
- }*/
-
// straight wall cell bounce back
- for ( std::size_t x = 2; x < dimX - 2; ++x ) {
- computeWallCell({x, 1}, { 0, 1});
- computeZouHeVelocityWallCell({x, dimY-2}, { 0,-1}, lidVelocityX);
+ for ( std::size_t x = 0; x < dimX; ++x ) {
+ computeWallCell({x, 0}, { 0, 1});
+ computeZouHeVelocityWallCell({x, dimY-1}, { 0,-1}, lidVelocityX);
+ //computeSimpleVelocityWallCell({x, dimY-1}, { 0,-1}, lidVelocityX);
}
- for ( std::size_t y = 2; y < dimY - 2; ++y ) {
- computeWallCell({1, y}, { 1, 0});
- computeWallCell({dimX-2, y}, {-1, 0});
+ for ( std::size_t y = 1; y < dimY-1; ++y ) {
+ computeWallCell({0, y}, { 1, 0});
+ computeWallCell({dimX-1, y}, {-1, 0});
}
// edge wall cell bounce back
- computeWallCell({1, 1 }, { 1, 1});
- computeWallCell({1, dimY-2}, { 1,-1});
- computeWallCell({dimX-2, 1 }, {-1, 1});
- computeWallCell({dimX-2, dimY-2}, {-1,-1});
+ //computeWallCell({0, dimY-1}, { 1,-1});
+ //computeWallCell({dimX-1, dimY-1}, {-1,-1});
+ //computeWallCell({0, 0 }, { 1, 1});
+ //computeWallCell({dimX-1, 0 }, {-1, 1});
// update density, velocity field
for ( std::size_t x = 0; x < dimX; ++x ) {
@@ -229,9 +219,8 @@ void computeLbmStep(std::size_t t) {
void init() {
for ( std::size_t x = 0; x < dimX; ++x ) {
for ( std::size_t y = 0; y < dimY; ++y ) {
- density[x][y] = 1.0;
velocity[x][y] = { 0.0, 0.0 };
- force[x][y] = { 0.0, 0.0 };
+ density[x][y] = 1.0;
pop.curr(x,y).equilibrize(density[x][y], velocity[x][y]);
pop.prev(x,y).equilibrize(density[x][y], velocity[x][y]);
@@ -247,33 +236,33 @@ void writeCurrentStateAsVTK(int time) {
fout << "lbm_output\n";
fout << "ASCII\n";
fout << "DATASET RECTILINEAR_GRID\n";
- fout << "DIMENSIONS " << dimX - 2 << " " << dimY - 2 << " 1" << "\n";
+ fout << "DIMENSIONS " << dimX << " " << dimY << " 1" << "\n";
- fout << "X_COORDINATES " << dimX - 2 << " float\n";
- for( std::size_t x = 1; x < dimX - 1; ++x ) {
+ fout << "X_COORDINATES " << dimX << " float\n";
+ for( std::size_t x = 0; x < dimX; ++x ) {
fout << x << " ";
}
- fout << "\nY_COORDINATES " << dimY - 2 << " float\n";
- for( std::size_t y = 1; y < dimY - 1; ++y ) {
+ fout << "\nY_COORDINATES " << dimY << " float\n";
+ for( std::size_t y = 0; y < dimY; ++y ) {
fout << y << " ";
}
fout << "\nZ_COORDINATES " << 1 << " float\n";
fout << 0 << "\n";
- fout << "POINT_DATA " << (dimX - 2) * (dimY - 2) << "\n";
+ fout << "POINT_DATA " << dimX * dimY << "\n";
fout << "VECTORS velocity float\n";
- for ( std::size_t y = 1; y < dimY - 1; ++y ) {
- for ( std::size_t x = 1; x < dimX - 1; ++x ) {
+ for ( std::size_t y = 0; y < dimY; ++y ) {
+ for ( std::size_t x = 0; x < dimX; ++x ) {
fout << velocity[x][y][0] << " " << velocity[x][y][1] << " 0\n";
}
}
fout << "SCALARS density float 1\n";
fout << "LOOKUP_TABLE default\n";
- for ( std::size_t y = 1; y < dimY - 1; ++y ) {
- for ( std::size_t x = 1; x < dimX - 1; ++x ) {
+ for ( std::size_t y = 0; y < dimY; ++y ) {
+ for ( std::size_t x = 0; x < dimX; ++x ) {
fout << density[x][y] << "\n";
}
}
@@ -284,13 +273,13 @@ void writeCurrentStateAsVTK(int time) {
int main() {
init();
- std::cout << "u: " << lidVelocityX << std::endl;
std::cout << "tau: " << tau << std::endl;
+ std::cout << "u: " << lidVelocityX << std::endl;
- for ( std::size_t t = 0; t < 20000; ++t ) {
+ for ( std::size_t t = 0; t <= 40000; ++t ) {
computeLbmStep(t);
- if ( t % 1000 == 0 ) {
+ if ( t < 100 || t % 1000 == 0 ) {
std::cout << "Writing " << t << "." << std::endl;
writeCurrentStateAsVTK(t);
}