Add inflow to channel

…it seems that as far as pure visuals are concerned simply not doing
anything for the outflow works out ok.

Actual outflow condition to complement this primitive Dirichlet
inflow remains to be implemented.
to be implemented.

1 files changed, 22 insertions, 33 deletions

diff --git a/channel.cc b/channel.cc
index 837b355..4307dfd 100644
--- a/channel.cc
+++ b/channel.cc
@@ -9,18 +9,19 @@
 constexpr std::size_t dimX = 500;
 constexpr std::size_t dimY = 40;
 
-constexpr double uWall    = 0.2;
-constexpr double reynolds = 500;
+constexpr double uInflow  = 0.1;
+constexpr double reynolds = 1000;
 
-constexpr double tau   = 3. * uWall * (dimX-1) / reynolds + 0.5;
+constexpr double tau   = 3. * uInflow * (dimX-1) / reynolds + 0.5;
 constexpr double omega = 1. / tau;
 
 DataCellBuffer pop(dimX, dimY);
 FluidBuffer fluid(dimX, dimY);
 
 std::vector<BoxObstacle> obstacles{
-	{300, 0,  320, 25},
-	{340, 15, 360, 39},
+	{100, 0,  120, 25},
+	{140, 15, 160, 39},
+	{300, 15, 340, 25},
 };
 
 void init() {
@@ -50,25 +51,7 @@ void computeLbmStep() {
 		}
 	}
 
-	// periodic boundary
-	for ( std::size_t y = 1; y < dimY-1; ++y ) {
-		pop.curr(0,y+1).get(1, 1) = pop.prev(dimX-1,y).get( 1, 1);
-		pop.curr(0,y  ).get(1, 0) = pop.prev(dimX-1,y).get( 1, 0);
-		pop.curr(0,y-1).get(1,-1) = pop.prev(dimX-1,y).get( 1,-1);
-	}
-	for ( std::size_t y = 1; y < dimY-1; ++y ) {
-		pop.curr(dimX-1,y+1).get(-1, 1) = pop.prev(0,y).get(-1, 1);
-		pop.curr(dimX-1,y  ).get(-1, 0) = pop.prev(0,y).get(-1, 0);
-		pop.curr(dimX-1,y-1).get(-1,-1) = pop.prev(0,y).get(-1,-1);
-	}
-
-	// straight wall cell bounce back
-	for ( std::size_t x = 0; x < 100; ++x ) {
-		computeZouHeVelocityWallCell(pop, {x, 0     }, { 0, 1}, uWall);
-		computeZouHeVelocityWallCell(pop, {x, dimY-1}, { 0,-1}, uWall);
-	}
-
-	for ( std::size_t x = 100; x < dimX-1; ++x ) {
+	for ( std::size_t x = 0; x < dimX; ++x ) {
 		computeWallCell(pop, {x, 0     }, { 0, 1});
 		computeWallCell(pop, {x, dimY-1}, { 0,-1});
 	}
@@ -80,9 +63,15 @@ void computeLbmStep() {
 
 	for ( std::size_t x = 0; x < dimX; ++x ) {
 		for ( std::size_t y = 0; y < dimY; ++y ) {
-			Cell& cell = static_cast<Cell&>(pop.curr(x,y));
-			fluid.density(x,y)  = cell.sum();
-			fluid.velocity(x,y) = cell.velocity(fluid.density(x,y));
+			if ( x == 0 ) {
+				// inflow
+				fluid.density(x,y)  = 1.0;
+				fluid.velocity(x,y) = { uInflow, 0 };
+			} else {
+				Cell& cell = static_cast<Cell&>(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) {
 				return !o.isInside(x, y);
@@ -96,15 +85,15 @@ void computeLbmStep() {
 int main() {
 	init();
 
-	std::cout << "Re:    " << reynolds << std::endl;
-	std::cout << "uWall: " << uWall    << std::endl;
-	std::cout << "tau:   " << tau      << std::endl;
-	std::cout << "omega: " << omega    << std::endl;
+	std::cout << "Re:      " << reynolds << std::endl;
+	std::cout << "uInflow: " << uInflow  << std::endl;
+	std::cout << "tau:     " << tau      << std::endl;
+	std::cout << "omega:   " << omega    << std::endl;
 
-	for ( std::size_t t = 0; t <= 10000; ++t ) {
+	for ( std::size_t t = 0; t <= 5000; ++t ) {
 		computeLbmStep();
 
-		if ( t % 100 == 0 ) {
+		if ( t % 50 == 0 ) {
 			std::cout << ".";
 			std::cout.flush();
 			fluid.writeAsVTK("result/data_t" + std::to_string(t) + ".vtk");