From 65aca58ba4167a8e4dfeb7274a5a43fb7343b2f8 Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Thu, 11 Oct 2018 22:34:42 +0200
Subject: Setup basic periodic channel example

---
 CMakeLists.txt                      |  12 +++-
 channel.cc                          | 115 ++++++++++++++++++++++++++++++++++++
 lid_driven_cavity.cc                |   3 -
 lid_driven_cavity_with_obstacles.cc |   1 -
 4 files changed, 126 insertions(+), 5 deletions(-)
 create mode 100644 channel.cc

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4248a32..bf8dc7b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -6,6 +6,7 @@ add_library(
 		src/lbm.cc
 		src/fluid_buffer.cc
 		src/boundary_conditions.cc
+		src/box_obstacle.cc
 )
 
 target_compile_features(
@@ -39,10 +40,19 @@ target_link_libraries(
 add_executable(
 	lid_driven_cavity_with_obstacles
 		lid_driven_cavity_with_obstacles.cc
-		src/box_obstacle.cc
 )
 
 target_link_libraries(
 	lid_driven_cavity_with_obstacles
 		boltzbub
 )
+
+add_executable(
+	channel
+		channel.cc
+)
+
+target_link_libraries(
+	channel
+		boltzbub
+)
diff --git a/channel.cc b/channel.cc
new file mode 100644
index 0000000..837b355
--- /dev/null
+++ b/channel.cc
@@ -0,0 +1,115 @@
+#include <iostream>
+#include <vector>
+#include <algorithm>
+
+#include "lbm.h"
+#include "boundary_conditions.h"
+#include "box_obstacle.h"
+
+constexpr std::size_t dimX = 500;
+constexpr std::size_t dimY = 40;
+
+constexpr double uWall    = 0.2;
+constexpr double reynolds = 500;
+
+constexpr double tau   = 3. * uWall * (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},
+};
+
+void init() {
+	for ( std::size_t x = 0; x < dimX; ++x ) {
+		for ( std::size_t y = 0; y < dimY; ++y ) {
+			fluid.density(x,y)  = 1.0;
+			fluid.velocity(x,y) = { 0.0, 0.0 };
+
+			static_cast<Cell&>(pop.curr(x,y)).equilibrize(
+				fluid.density(x,y), fluid.velocity(x,y));
+			static_cast<Cell&>(pop.prev(x,y)).equilibrize(
+				fluid.density(x,y), fluid.velocity(x,y));
+		}
+	}
+}
+
+void computeLbmStep() {
+	pop.swap();
+
+	for ( std::size_t x = 0; x < dimX; ++x ) {
+		for ( std::size_t y = 0; y < dimY; ++y ) {
+			if ( std::all_of(obstacles.cbegin(), obstacles.cend(), [x, y](const auto& o) {
+				return !o.isInside(x, y);
+			}) ) {
+				streamFluidCell(pop, x, y);
+			}
+		}
+	}
+
+	// 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 ) {
+		computeWallCell(pop, {x, 0     }, { 0, 1});
+		computeWallCell(pop, {x, dimY-1}, { 0,-1});
+	}
+
+	// obstacles
+	for ( const auto& box : obstacles ) {
+		box.applyBoundary(pop);
+	}
+
+	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 ( std::all_of(obstacles.cbegin(), obstacles.cend(), [x, y](const auto& o) {
+				return !o.isInside(x, y);
+			}) ) {
+				collideFluidCell(omega, pop, fluid, x, y);
+			}
+		}
+	}
+}
+
+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;
+
+	for ( std::size_t t = 0; t <= 10000; ++t ) {
+		computeLbmStep();
+
+		if ( t % 100 == 0 ) {
+			std::cout << ".";
+			std::cout.flush();
+			fluid.writeAsVTK("result/data_t" + std::to_string(t) + ".vtk");
+		}
+	}
+
+	std::cout << std::endl;
+}
diff --git a/lid_driven_cavity.cc b/lid_driven_cavity.cc
index ca4b030..79f330b 100644
--- a/lid_driven_cavity.cc
+++ b/lid_driven_cavity.cc
@@ -1,7 +1,4 @@
 #include <iostream>
-#include <vector>
-#include <string>
-#include <algorithm>
 
 #include "lbm.h"
 #include "boundary_conditions.h"
diff --git a/lid_driven_cavity_with_obstacles.cc b/lid_driven_cavity_with_obstacles.cc
index 092814e..55f519e 100644
--- a/lid_driven_cavity_with_obstacles.cc
+++ b/lid_driven_cavity_with_obstacles.cc
@@ -1,6 +1,5 @@
 #include <iostream>
 #include <vector>
-#include <string>
 #include <algorithm>
 
 #include "lbm.h"
-- 
cgit v1.2.3