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 --- 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 ++ 7 files changed, 331 insertions(+) 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 (limited to 'src') 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