From acdba1d0ba0de53df326956352ff09bea2b80c6a Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Sun, 9 Sep 2018 12:08:26 +0200 Subject: Split into compilation units --- CMakeLists.txt | 4 +++ src/boundary_conditions.cc | 45 ++++++++++++++++++++++++ src/boundary_conditions.h | 44 ++--------------------- src/box_obstacle.cc | 28 +++++++++++++++ src/box_obstacle.h | 26 ++------------ src/data_cell_buffer.h | 82 +++++++++++++++++++++--------------------- src/fluid_buffer.cc | 75 +++++++++++++++++++++++++++++++++++++++ src/fluid_buffer.h | 88 ++++++++-------------------------------------- src/lbm.cc | 68 +++++++++++++++++++++++++++++++++++ src/lbm.h | 65 +++------------------------------- 10 files changed, 285 insertions(+), 240 deletions(-) create mode 100644 src/boundary_conditions.cc create mode 100644 src/box_obstacle.cc create mode 100644 src/fluid_buffer.cc create mode 100644 src/lbm.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index a5b185d..8818966 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,4 +13,8 @@ include_directories( add_executable( boltz boltz.cc + src/lbm.cc + src/fluid_buffer.cc + src/boundary_conditions.cc + src/box_obstacle.cc ) diff --git a/src/boundary_conditions.cc b/src/boundary_conditions.cc new file mode 100644 index 0000000..3488463 --- /dev/null +++ b/src/boundary_conditions.cc @@ -0,0 +1,45 @@ +#include "boundary_conditions.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/boundary_conditions.h b/src/boundary_conditions.h index c825c32..6b7de19 100644 --- a/src/boundary_conditions.h +++ b/src/boundary_conditions.h @@ -3,46 +3,6 @@ #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); -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 ); -} +void computeZouHeVelocityWallCell(DataCellBuffer& pop, Vector cell, Vector normal, double vX); diff --git a/src/box_obstacle.cc b/src/box_obstacle.cc new file mode 100644 index 0000000..9cb6ad0 --- /dev/null +++ b/src/box_obstacle.cc @@ -0,0 +1,28 @@ +#include "box_obstacle.h" + +#include "boundary_conditions.h" + +BoxObstacle::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 BoxObstacle::isInside(std::size_t x, std::size_t y) const { + return x > lower_x_ + && x < upper_x_ + && y > lower_y_ + && y < upper_y_; +} + +void BoxObstacle::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/box_obstacle.h b/src/box_obstacle.h index 46660a7..149e8e1 100644 --- a/src/box_obstacle.h +++ b/src/box_obstacle.h @@ -1,7 +1,6 @@ #pragma once #include "data_cell_buffer.h" -#include "boundary_conditions.h" struct BoxObstacle { const std::size_t lower_x_; @@ -9,28 +8,9 @@ struct BoxObstacle { 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) { } + BoxObstacle(std::size_t lX, std::size_t lY, std::size_t uX, std::size_t 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_; - } + bool isInside(std::size_t x, std::size_t y) const; - 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}); - } + void applyBoundary(DataCellBuffer& pop) const; }; diff --git a/src/data_cell_buffer.h b/src/data_cell_buffer.h index 06c2644..033fb24 100644 --- a/src/data_cell_buffer.h +++ b/src/data_cell_buffer.h @@ -5,46 +5,46 @@ #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]); - } +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.cc b/src/fluid_buffer.cc new file mode 100644 index 0000000..e73c818 --- /dev/null +++ b/src/fluid_buffer.cc @@ -0,0 +1,75 @@ +#include "fluid_buffer.h" + +#include + +FluidBuffer::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 FluidBuffer::dimX() const { + return dim_x_; +} + +std::size_t FluidBuffer::dimY() const { + return dim_y_; +} + +Density& FluidBuffer::density(std::size_t x, std::size_t y) { + return density_[y*dim_x_ + x]; +} + +Density& FluidBuffer::density(Vector pos) { + return density(pos[0], pos[1]); +} + +Velocity& FluidBuffer::velocity(std::size_t x, std::size_t y) { + return velocity_[y*dim_x_ + x]; +} + +Velocity& FluidBuffer::velocity(Vector pos) { + return velocity(pos[0], pos[1]); +} + +void FluidBuffer::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/fluid_buffer.h b/src/fluid_buffer.h index a6c41ba..1d802c0 100644 --- a/src/fluid_buffer.h +++ b/src/fluid_buffer.h @@ -1,6 +1,5 @@ #pragma once -#include #include #include "vector.h" @@ -9,83 +8,24 @@ using Velocity = Vector; using Density = double; class FluidBuffer { - private: - const std::size_t dim_x_; - const std::size_t dim_y_; +private: + const std::size_t dim_x_; + const std::size_t dim_y_; - std::unique_ptr density_; - std::unique_ptr velocity_; + 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]) { } +public: + FluidBuffer(std::size_t dimX, std::size_t dimY); - std::size_t dimX() const { - return dim_x_; - } + std::size_t dimX() const; + std::size_t dimY() const; - std::size_t dimY() const { - return dim_y_; - } + Density& density(std::size_t x, std::size_t y); + Density& density(Vector pos); - Density& density(std::size_t x, std::size_t y) { - return density_[y*dim_x_ + x]; - } + Velocity& velocity(std::size_t x, std::size_t y); + Velocity& velocity(Vector pos); - 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(); - } + void writeAsVTK(int time); }; diff --git a/src/lbm.cc b/src/lbm.cc new file mode 100644 index 0000000..ca84cee --- /dev/null +++ b/src/lbm.cc @@ -0,0 +1,68 @@ +#include "lbm.h" + +void Cell::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 Cell::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 Cell::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/lbm.h b/src/lbm.h index 685fa08..96348c1 100644 --- a/src/lbm.h +++ b/src/lbm.h @@ -11,71 +11,16 @@ constexpr DataCell weight{ }; 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())); - } - } - } + void equilibrize(Density d, Velocity v); - 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); - } + double sum() const; - 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) - }; - } + Velocity velocity(Density d) const; }; -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 streamFluidCell(DataCellBuffer& pop, std::size_t x, std::size_t y); 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)); - } - } - } - } -} - + std::size_t x, std::size_t y); -- cgit v1.2.3