From 1a30c8427f7711064462ee436de614f90dbad002 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Sat, 31 Mar 2018 18:09:09 +0200 Subject: Fix and improve wall handling Use three neighbor cells for bounceback instead of three directions of a single cell. Generalize neighbor cell determination using normal vector of wall. --- boltz.cc | 113 ++++++++++++++++++++++++++++++++++------------------------- src/vector.h | 22 ++++++------ 2 files changed, 78 insertions(+), 57 deletions(-) diff --git a/boltz.cc b/boltz.cc index 85d939a..1acebc2 100644 --- a/boltz.cc +++ b/boltz.cc @@ -18,8 +18,8 @@ struct DataCell { } }; -using Velocity = Vector; -using Force = Vector; +using Velocity = Vector; +using Force = Vector; using Density = double; constexpr DataCell weight{ @@ -81,7 +81,7 @@ class CellBuffer { constexpr std::size_t dimX = 128; constexpr std::size_t dimY = 128; -constexpr double tau = 0.6; +constexpr double tau = 0.8; CellBuffer pop(dimX, dimY); @@ -90,6 +90,8 @@ Velocity velocity[dimX][dimY]; Force force [dimX][dimY]; void computeFluidCell(std::size_t x, std::size_t y) { + //velocity[x][y] += tau / density[x][y] * force[x][y]; + // compute equilibrium Cell eq; eq.equilibrize(density[x][y], velocity[x][y]); @@ -102,14 +104,31 @@ void computeFluidCell(std::size_t x, std::size_t y) { } } -inline int clamp(int x) { - return x / -x; +std::pair, Vector> neighbors(Vector v) { + if ( v[0] == 0 ) { + return std::pair, Vector>{ + { -1, v[1] }, + { 1, v[1] } + }; + } else if ( v[1] == 0 ) { + return std::pair, Vector>{ + { v[0], -1 }, + { v[0], 1 } + }; + } else { + return std::pair, Vector>{ + { 0, v[1] }, + { v[0], 0 } + }; + } } -void computeWallCell(std::size_t x, std::size_t y, int normalX, int normalY) { - pop.curr(x,y).get(clamp(normalX+normalY),clamp(normalY-normalX)) = pop.curr(x-normalX,y-normalY).get(clamp(normalX-normalY),clamp(normalY+normalX)); - pop.curr(x,y).get(normalX ,normalY ) = pop.curr(x-normalX,y-normalY).get(-normalX,-normalY); - pop.curr(x,y).get(clamp(normalX-normalY),clamp(normalY+normalX)) = pop.curr(x-normalX,y-normalY).get(clamp(normalX+normalY),clamp(normalY-normalX)); +void computeWallCell(std::size_t x, std::size_t y, Vector normal) { + const auto [neighborA, neighborB] = neighbors(normal); + + pop.curr(x,y).get(neighborA[0], neighborA[1]) = pop.curr(x-neighborA[0], y-neighborA[1]).get(-neighborA[0], -neighborA[1]); + pop.curr(x,y).get(normal[0] , normal[1] ) = pop.curr(x-normal[0] , y-normal[1] ).get(-normal[0] , -normal[1] ); + pop.curr(x,y).get(neighborB[0], neighborB[1]) = pop.curr(x-neighborB[0], y-neighborB[1]).get(-neighborB[0], -neighborB[1]); } void computeLbmStep(std::size_t t) { @@ -117,7 +136,7 @@ void computeLbmStep(std::size_t t) { 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 >= 40 ) { + if ( x <= 20 || x >= 40 || y <= 20 || y >= 80 ) { computeFluidCell(x, y); } } @@ -126,35 +145,35 @@ void computeLbmStep(std::size_t t) { // obstacle { for ( std::size_t x = 21; x < 40; ++x ) { - computeWallCell(x, 40, 0, 1); - computeWallCell(x, 20, 0, -1); + computeWallCell(x, 80, { 0, 1}); + computeWallCell(x, 20, { 0,-1}); } - for ( std::size_t y = 21; y < 40; ++y ) { - computeWallCell(40, y, 1, 0); - computeWallCell(20, y, -1, 0); + for ( std::size_t y = 21; y < 80; ++y ) { + computeWallCell(40, y, { 1, 0}); + computeWallCell(20, y, {-1, 0}); } - computeWallCell(40,40, 1, 1); - computeWallCell(40,20, 1,-1); - computeWallCell(20,40,-1, 1); - computeWallCell(20,20,-1,-1); + 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); - computeWallCell(x, dimY-2, 0, -1); + computeWallCell(x, 1, { 0, 1}); + computeWallCell(x, dimY-2,{ 0,-1}); } for ( std::size_t y = 2; y < dimY - 2; ++y ) { - computeWallCell(1,y,1,0); - computeWallCell(dimX-2,y,-1,0); + computeWallCell(1, y, { 1, 0}); + computeWallCell(dimX-2, 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(1, 1, { 1, 1}); + computeWallCell(1, dimY-2, { 1,-1}); + computeWallCell(dimX-2, 1, {-1, 1}); + computeWallCell(dimX-2, dimY-2, {-1,-1}); // update density, velocity field for ( std::size_t x = 0; x < dimX; ++x ) { @@ -165,6 +184,25 @@ 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 }; + + pop.curr(x,y).equilibrize(density[x][y], velocity[x][y]); + pop.prev(x,y).equilibrize(density[x][y], velocity[x][y]); + } + } + + for ( std::size_t y = 50; y < dimY-50; ++y ) { + for ( std::size_t x = 50; x < dimX-50; ++x ) { + density[x][y] = 2.0; + } + } +} + void writeCurrentStateAsVTK(int time) { std::ofstream fout; fout.open(("result/data_t" + std::to_string(time) + ".vtk").c_str()); @@ -207,29 +245,10 @@ void writeCurrentStateAsVTK(int time) { fout.close(); } -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 }; - - pop.curr(x,y).equilibrize(density[x][y], velocity[x][y]); - pop.prev(x,y).equilibrize(density[x][y], velocity[x][y]); - } - } - - for ( std::size_t y = 55; y < dimY-55; ++y ) { - for ( std::size_t x = 75; x < dimX-35; ++x ) { - density[x][y] = 0.6; - } - } -} - int main() { init(); - for ( std::size_t t = 0; t < 500; ++t ) { + for ( std::size_t t = 0; t < 1000; ++t ) { computeLbmStep(t); if ( t % 2 == 0 ) { diff --git a/src/vector.h b/src/vector.h index 55c7663..e0dc239 100644 --- a/src/vector.h +++ b/src/vector.h @@ -4,41 +4,43 @@ inline double sq(double x) noexcept { return x * x; } +template struct Vector { - double data[2]; + T data[2]; - double comp(int x, int y) const { + T comp(int x, int y) const { return x*data[0] + y*data[1]; } - double norm() const { + T norm() const { return std::sqrt(sq(data[0]) + sq(data[1])); } - double& operator[](std::size_t i) { + T& operator[](std::size_t i) { return data[i]; } - double operator[](std::size_t i) const { + T operator[](std::size_t i) const { return data[i]; } - Vector operator*(double scalar) const { - return Vector{ + Vector operator*(T scalar) const { + return Vector{ data[0] * scalar, data[1] * scalar }; } - Vector& operator+=(const Vector& rhs) { + Vector& operator+=(const Vector& rhs) { data[0] += rhs[0]; data[1] += rhs[1]; return *this; } }; -Vector operator*(double scalar, const Vector& vec) { - return Vector{ +template +Vector operator*(T scalar, const Vector& vec) { + return Vector{ vec[0] * scalar, vec[1] * scalar }; -- cgit v1.2.3