From c4771584ae3b0b4c34fbef20f13c6fc7d80cd046 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Thu, 2 Aug 2018 21:07:30 +0200 Subject: Bring in some more structure --- boltz.cc | 322 ++++++---------------------------------------- src/boundary_conditions.h | 48 +++++++ src/box_obstacle.h | 36 ++++++ src/data_cell.h | 23 ++++ src/data_cell_buffer.h | 50 +++++++ src/fluid_buffer.h | 91 +++++++++++++ src/lbm.h | 81 ++++++++++++ src/vector.h | 2 + 8 files changed, 368 insertions(+), 285 deletions(-) create mode 100644 src/boundary_conditions.h create mode 100644 src/box_obstacle.h create mode 100644 src/data_cell.h create mode 100644 src/data_cell_buffer.h create mode 100644 src/fluid_buffer.h create mode 100644 src/lbm.h diff --git a/boltz.cc b/boltz.cc index e36e337..6f8788e 100644 --- a/boltz.cc +++ b/boltz.cc @@ -1,117 +1,10 @@ #include -#include -#include -#include -#include #include #include -#include "vector.h" - -struct DataCell { - double data[3][3]; - - inline double& get(int x, int y) { - return data[1+x][1-y]; - } - - inline double& get(Vector v) { - return get(v[0], v[1]); - } - - inline double get(int x, int y) const { - return data[1+x][1-y]; - } - - inline double get(Vector v) const { - return get(v[0], v[1]); - } -}; - -using Velocity = Vector; -using Density = double; - -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) - }; - } -}; - -class CellBuffer { - private: - const std::size_t dim_x_; - const std::size_t dim_y_; - - std::unique_ptr curr_; - std::unique_ptr 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& curr(Vector pos) { - return curr(pos[0], pos[1]); - } - - inline Cell& prev(std::size_t x, std::size_t y) { - return prev_[y*dim_x_ + x]; - } - - inline Cell& prev(Vector pos) { - return prev(pos[0], pos[1]); - } -}; - -std::pair, Vector> neighbors(Vector v) { - if ( v[0] == 0 ) { - return { - { -1, v[1] }, - { 1, v[1] } - }; - } else if ( v[1] == 0 ) { - return { - { v[0], -1 }, - { v[0], 1 } - }; - } else { - return { - { 0, v[1] }, - { v[0], 0 } - }; - } -} +#include "lbm.h" +#include "boundary_conditions.h" +#include "box_obstacle.h" constexpr std::size_t dimX = 120; constexpr std::size_t dimY = dimX; @@ -122,110 +15,8 @@ constexpr double reynolds = 1000; constexpr double tau = 3 * uLid * (dimX-1) / reynolds + 0.5; constexpr double omega = 1. / tau; -CellBuffer pop(dimX, dimY); - -Density density [dimX][dimY]; -Velocity velocity[dimX][dimY]; - -void streamFluidCell(std::size_t x, std::size_t y) { - if ( x != 0 && x != dimX - 1 && y != 0 && y != 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+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); - } - } - } - } -} - -void collideFluidCell(std::size_t x, std::size_t y) { - // compute equilibrium - Cell eq; - eq.equilibrize(density[x][y], velocity[x][y]); - - // collide (BGK, relax towards equilibrium) - 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,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+i >= 0 && x+i <= dimX - 1 && y+j >= 0 && y+j <= 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)); - } - } - } - } -} - -void computeWallCell(Vector cell, Vector normal) { - const auto [neighborA, neighborB] = neighbors(normal); - - 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 computeZouHeVelocityWallCell(Vector cell, Vector normal, double vX) { - const auto [neighborA, neighborB] = neighbors(normal); - - 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).get(-neighborA) + - pop.curr(cell).get(-normal ) + - pop.curr(cell).get(-neighborB) - ); - - 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); -} - -struct BoxObstacle { - const std::size_t lower_x_; - const std::size_t lower_y_; - const std::size_t upper_x_; - const std::size_t upper_y_; - - BoxObstacle(std::size_t lX, std::size_t lY, std::size_t uX, std::size_t uY): - lower_x_(lX), lower_y_(lY), upper_x_(uX), upper_y_(uY) { } - - bool isInside(std::size_t x, std::size_t y) const { - return x > lower_x_ - && x < upper_x_ - && y > lower_y_ - && y < upper_y_; - } - - void applyBoundary() const { - for ( std::size_t x = lower_x_+1; x < upper_x_; ++x ) { - computeWallCell({x, lower_y_}, { 0,-1}); - computeWallCell({x, upper_y_}, { 0, 1}); - } - for ( std::size_t y = lower_y_+1; y < upper_y_; ++y ) { - computeWallCell({lower_x_, y}, {-1, 0}); - computeWallCell({upper_x_, y}, { 1, 0}); - } - computeWallCell({lower_x_, lower_y_}, {-1,-1}); - computeWallCell({upper_x_, lower_y_}, { 1,-1}); - computeWallCell({upper_x_, upper_y_}, { 1, 1}); - computeWallCell({lower_x_, upper_y_}, {-1, 1}); - } -}; +DataCellBuffer pop(dimX, dimY); +FluidBuffer fluid(dimX, dimY); std::vector obstacles{ {20, 20, 40, 40}, @@ -239,108 +30,69 @@ std::vector obstacles{ {80, 80, 100, 100}, }; -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 ) { + fluid.density(x,y) = 1.0; + fluid.velocity(x,y) = { 0.0, 0.0 }; + + static_cast(pop.curr(x,y)).equilibrize( + fluid.density(x,y), fluid.velocity(x,y)); + static_cast(pop.prev(x,y)).equilibrize( + fluid.density(x,y), fluid.velocity(x,y)); + } + } +} + +void computeLbmStep() { pop.swap(); for ( std::size_t x = 0; x < dimX; ++x ) { for ( std::size_t y = 0; y < dimY; ++y ) { - if ( std::all_of(obstacles.cbegin(), obstacles.cend(), [x, y](const auto& o) -> bool { + if ( std::all_of(obstacles.cbegin(), obstacles.cend(), [x, y](const auto& o) { return !o.isInside(x, y); }) ) { - streamFluidCell(x, y); + streamFluidCell(pop, x, y); } } } // straight wall cell bounce back for ( std::size_t x = 0; x < dimX; ++x ) { - computeZouHeVelocityWallCell({x, dimY-1}, { 0,-1}, uLid); + computeZouHeVelocityWallCell(pop, {x, dimY-1}, { 0,-1}, uLid); } for ( std::size_t x = 1; x < dimX-1; ++x ) { - computeWallCell({x, 0}, { 0, 1}); + computeWallCell(pop, {x, 0}, { 0, 1}); } for ( std::size_t y = 1; y < dimY-1; ++y ) { - computeWallCell({0, y}, { 1, 0}); - computeWallCell({dimX-1, y}, {-1, 0}); + computeWallCell(pop, {0, y}, { 1, 0}); + computeWallCell(pop, {dimX-1, y}, {-1, 0}); } // edge wall cell bounce back - computeWallCell({0, 0 }, { 1, 1}); - computeWallCell({dimX-1, 0 }, {-1, 1}); + computeWallCell(pop, {0, 0 }, { 1, 1}); + computeWallCell(pop, {dimX-1, 0 }, {-1, 1}); // obstacles - for ( auto& o : obstacles ) { - o.applyBoundary(); + for ( const auto& box : obstacles ) { + box.applyBoundary(pop); } for ( std::size_t x = 0; x < dimX; ++x ) { for ( std::size_t y = 0; y < dimY; ++y ) { - density[x][y] = pop.curr(x,y).sum(); - velocity[x][y] = pop.curr(x,y).velocity(density[x][y]); + Cell& cell = static_cast(pop.curr(x,y)); + fluid.density(x,y) = cell.sum(); + fluid.velocity(x,y) = cell.velocity(fluid.density(x,y)); - if ( std::all_of(obstacles.cbegin(), obstacles.cend(), [x, y](const auto& o) -> bool { + if ( std::all_of(obstacles.cbegin(), obstacles.cend(), [x, y](const auto& o) { return !o.isInside(x, y); }) ) { - collideFluidCell(x, y); + collideFluidCell(omega, pop, fluid, x, y); } } } } -void init() { - for ( std::size_t x = 0; x < dimX; ++x ) { - for ( std::size_t y = 0; y < dimY; ++y ) { - velocity[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]); - } - } -} - -void writeCurrentStateAsVTK(int time) { - std::ofstream fout; - fout.open(("result/data_t" + std::to_string(time) + ".vtk").c_str()); - - fout << "# vtk DataFile Version 3.0\n"; - fout << "lbm_output\n"; - fout << "ASCII\n"; - fout << "DATASET RECTILINEAR_GRID\n"; - fout << "DIMENSIONS " << dimX << " " << dimY << " 1" << "\n"; - - fout << "X_COORDINATES " << dimX << " float\n"; - for( std::size_t x = 0; x < dimX; ++x ) { - fout << x << " "; - } - - 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 * dimY << "\n"; - - fout << "VECTORS velocity float\n"; - 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 = 0; y < dimY; ++y ) { - for ( std::size_t x = 0; x < dimX; ++x ) { - fout << density[x][y] << "\n"; - } - } - - fout.close(); -} - int main() { init(); @@ -349,12 +101,12 @@ int main() { std::cout << "tau: " << tau << std::endl; for ( std::size_t t = 0; t <= 6000; ++t ) { - computeLbmStep(t); + computeLbmStep(); if ( t % 100 == 0 ) { std::cout << "."; std::cout.flush(); - writeCurrentStateAsVTK(t); + fluid.writeAsVTK(t); } } diff --git a/src/boundary_conditions.h b/src/boundary_conditions.h new file mode 100644 index 0000000..c825c32 --- /dev/null +++ b/src/boundary_conditions.h @@ -0,0 +1,48 @@ +#pragma once + +#include "vector.h" +#include "data_cell_buffer.h" + +std::pair, Vector> neighbors(Vector v) { + if ( v[0] == 0 ) { + return { + { -1, v[1] }, + { 1, v[1] } + }; + } else if ( v[1] == 0 ) { + return { + { v[0], -1 }, + { v[0], 1 } + }; + } else { + return { + { 0, v[1] }, + { v[0], 0 } + }; + } +} + +void computeWallCell(DataCellBuffer& pop, Vector cell, Vector normal) { + const auto [neighborA, neighborB] = neighbors(normal); + + 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 computeZouHeVelocityWallCell(DataCellBuffer& pop, Vector cell, Vector normal, double vX) { + const auto [neighborA, neighborB] = neighbors(normal); + + 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).get(-neighborA) + + pop.curr(cell).get(-normal ) + + pop.curr(cell).get(-neighborB) + ); + + 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 ); +} diff --git a/src/box_obstacle.h b/src/box_obstacle.h new file mode 100644 index 0000000..46660a7 --- /dev/null +++ b/src/box_obstacle.h @@ -0,0 +1,36 @@ +#pragma once + +#include "data_cell_buffer.h" +#include "boundary_conditions.h" + +struct BoxObstacle { + const std::size_t lower_x_; + const std::size_t lower_y_; + const std::size_t upper_x_; + const std::size_t upper_y_; + + BoxObstacle(std::size_t lX, std::size_t lY, std::size_t uX, std::size_t uY): + lower_x_(lX), lower_y_(lY), upper_x_(uX), upper_y_(uY) { } + + bool isInside(std::size_t x, std::size_t y) const { + return x > lower_x_ + && x < upper_x_ + && y > lower_y_ + && y < upper_y_; + } + + void applyBoundary(DataCellBuffer& pop) const { + for ( std::size_t x = lower_x_+1; x < upper_x_; ++x ) { + computeWallCell(pop, {x, lower_y_}, { 0,-1}); + computeWallCell(pop, {x, upper_y_}, { 0, 1}); + } + for ( std::size_t y = lower_y_+1; y < upper_y_; ++y ) { + computeWallCell(pop, {lower_x_, y}, {-1, 0}); + computeWallCell(pop, {upper_x_, y}, { 1, 0}); + } + computeWallCell(pop, {lower_x_, lower_y_}, {-1,-1}); + computeWallCell(pop, {upper_x_, lower_y_}, { 1,-1}); + computeWallCell(pop, {upper_x_, upper_y_}, { 1, 1}); + computeWallCell(pop, {lower_x_, upper_y_}, {-1, 1}); + } +}; diff --git a/src/data_cell.h b/src/data_cell.h new file mode 100644 index 0000000..fbb5dfb --- /dev/null +++ b/src/data_cell.h @@ -0,0 +1,23 @@ +#pragma once + +#include "vector.h" + +struct DataCell { + double data[3][3]; + + inline double& get(int x, int y) { + return data[1+x][1-y]; + } + + inline double& get(Vector v) { + return get(v[0], v[1]); + } + + inline double get(int x, int y) const { + return data[1+x][1-y]; + } + + inline double get(Vector v) const { + return get(v[0], v[1]); + } +}; diff --git a/src/data_cell_buffer.h b/src/data_cell_buffer.h new file mode 100644 index 0000000..06c2644 --- /dev/null +++ b/src/data_cell_buffer.h @@ -0,0 +1,50 @@ +#pragma once + +#include + +#include "data_cell.h" + +class DataCellBuffer { + private: + const std::size_t dim_x_; + const std::size_t dim_y_; + + std::unique_ptr curr_; + std::unique_ptr prev_; + + public: + DataCellBuffer(std::size_t dimX, std::size_t dimY): + dim_x_(dimX), + dim_y_(dimY), + curr_(new DataCell[dimX*dimY]), + prev_(new DataCell[dimX*dimY]) { } + + std::size_t dimX() const { + return dim_x_; + } + + std::size_t dimY() const { + return dim_y_; + } + + void swap() { + curr_.swap(prev_); + } + + inline DataCell& curr(std::size_t x, std::size_t y) { + return curr_[y*dim_x_ + x]; + } + + inline DataCell& curr(Vector pos) { + return curr(pos[0], pos[1]); + } + + inline DataCell& prev(std::size_t x, std::size_t y) { + return prev_[y*dim_x_ + x]; + } + + inline DataCell& prev(Vector pos) { + return prev(pos[0], pos[1]); + } +}; + diff --git a/src/fluid_buffer.h b/src/fluid_buffer.h new file mode 100644 index 0000000..a6c41ba --- /dev/null +++ b/src/fluid_buffer.h @@ -0,0 +1,91 @@ +#pragma once + +#include +#include + +#include "vector.h" + +using Velocity = Vector; +using Density = double; + +class FluidBuffer { + private: + const std::size_t dim_x_; + const std::size_t dim_y_; + + std::unique_ptr density_; + std::unique_ptr velocity_; + + public: + FluidBuffer(std::size_t dimX, std::size_t dimY): + dim_x_(dimX), + dim_y_(dimY), + density_(new Density[dimX*dimY]), + velocity_(new Velocity[dimX*dimY]) { } + + std::size_t dimX() const { + return dim_x_; + } + + std::size_t dimY() const { + return dim_y_; + } + + Density& density(std::size_t x, std::size_t y) { + return density_[y*dim_x_ + x]; + } + + Density& curr(Vector pos) { + return density(pos[0], pos[1]); + } + + Velocity& velocity(std::size_t x, std::size_t y) { + return velocity_[y*dim_x_ + x]; + } + + Velocity& velocity(Vector pos) { + return velocity(pos[0], pos[1]); + } + + void writeAsVTK(int time) { + std::ofstream fout; + fout.open(("result/data_t" + std::to_string(time) + ".vtk").c_str()); + + fout << "# vtk DataFile Version 3.0\n"; + fout << "lbm_output\n"; + fout << "ASCII\n"; + fout << "DATASET RECTILINEAR_GRID\n"; + fout << "DIMENSIONS " << dimX() << " " << dimY() << " 1" << "\n"; + + fout << "X_COORDINATES " << dimX() << " float\n"; + for( std::size_t x = 0; x < dimX(); ++x ) { + fout << x << " "; + } + + 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() * dimY() << "\n"; + + fout << "VECTORS velocity float\n"; + 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 = 0; y < dimY(); ++y ) { + for ( std::size_t x = 0; x < dimX(); ++x ) { + fout << density(x,y) << "\n"; + } + } + + fout.close(); + } +}; diff --git a/src/lbm.h b/src/lbm.h new file mode 100644 index 0000000..685fa08 --- /dev/null +++ b/src/lbm.h @@ -0,0 +1,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)); + } + } + } + } +} + diff --git a/src/vector.h b/src/vector.h index 3a3ed5c..5e105c4 100644 --- a/src/vector.h +++ b/src/vector.h @@ -1,3 +1,5 @@ +#pragma once + #include inline double sq(double x) noexcept { -- cgit v1.2.3