From 9ca651ff422e7da7a2c1ffdfec6954c9a478b553 Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Fri, 12 Oct 2018 22:20:09 +0200
Subject: 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.
---
 channel.cc | 55 ++++++++++++++++++++++---------------------------------
 1 file 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");
-- 
cgit v1.2.3