Implement moving wall / velocity Dirichlet boundary condition

Usable as both the inflow condition of the channel example and the top
wall of a lid driven cavity.

6 files changed, 48 insertions, 16 deletions

diff --git a/README.md b/README.md
index 385627b..5ab879a 100644
--- a/README.md
+++ b/README.md
@@ -10,7 +10,7 @@
 | -                                     | -                                                                                                              |
 | `lid_driven_cavity.cc`                | Lid driven cavity using Zou/He boundary conditions for the top wall and simple bounce back for all other walls |
 | `lid_driven_cavity_with_obstacles.cc` | Same as `lid_driven_cavity.cc` but includes a grid of boxes to make things more interesting                    |
-| `channel.cc`                          | Channel flow some obstacles and Dirichlet inflow                                                               |
+| `channel.cc`                          | Channel flow with some obstacles and Dirichlet inflow condition                                                |
 
 ## Build
 
diff --git a/channel.cc b/channel.cc
index 3e7641b..b1a6166 100644
--- a/channel.cc
+++ b/channel.cc
@@ -9,8 +9,8 @@
 constexpr std::size_t dimX = 500;
 constexpr std::size_t dimY = 40;
 
-constexpr double uInflow  = 0.1;
-constexpr double reynolds = 1000;
+constexpr double uInflow  = 0.01;
+constexpr double reynolds = 100;
 
 constexpr double tau   = 3. * uInflow * (dimX-1) / reynolds + 0.5;
 constexpr double omega = 1. / tau;
@@ -56,6 +56,10 @@ void computeLbmStep() {
 		computeWallCell(pop, {x, dimY-1}, { 0,-1});
 	}
 
+	for ( std::size_t y = 1; y < dimY-1; ++y ) {
+		computeMovingWallCell(pop, {0,y}, {1,0}, {uInflow,0});
+	}
+
 	// obstacles
 	for ( const auto& box : obstacles ) {
 		box.applyBoundary(pop);
@@ -63,16 +67,18 @@ void computeLbmStep() {
 
 	for ( std::size_t x = 0; x < dimX; ++x ) {
 		for ( std::size_t y = 0; y < dimY; ++y ) {
-			if ( x == 0 && y > 0 && y < dimY-1 ) {
-				// 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));
+			Cell& cell = static_cast<Cell&>(pop.curr(x,y));
+
+			// bulk density
+			fluid.density(x,y) = cell.sum();
+
+			// outflow
+			if ( x == dimX-1 && y > 0 && y < dimY-1 ) {
+				fluid.density(x,y) = 1.0;
 			}
 
+			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);
 			}) ) {
@@ -90,10 +96,10 @@ int main() {
 	std::cout << "tau:     " << tau      << std::endl;
 	std::cout << "omega:   " << omega    << std::endl;
 
-	for ( std::size_t t = 0; t <= 5000; ++t ) {
+	for ( std::size_t t = 0; t <= 10000; ++t ) {
 		computeLbmStep();
 
-		if ( t % 100 == 0 ) {
+		if ( t % 1000 == 0 ) {
 			std::cout << ".";
 			std::cout.flush();
 			fluid.writeAsVTK("result/data_t" + std::to_string(t) + ".vtk");
diff --git a/lid_driven_cavity.cc b/lid_driven_cavity.cc
index 79f330b..7a978ad 100644
--- a/lid_driven_cavity.cc
+++ b/lid_driven_cavity.cc
@@ -38,10 +38,12 @@ void computeLbmStep() {
 		}
 	}
 
-	// straight wall cell bounce back
+	// moving top wall
 	for ( std::size_t x = 0; x < dimX; ++x ) {
-		computeZouHeVelocityWallCell(pop, {x, dimY-1}, { 0,-1}, uLid);
+		computeMovingWallCell(pop, {x, dimY-1}, {0, -1}, {uLid, 0});
 	}
+
+	// straight wall cell bounce back
 	for ( std::size_t x = 1; x < dimX-1; ++x ) {
 		computeWallCell(pop, {x, 0}, { 0, 1});
 	}
@@ -76,7 +78,7 @@ int main() {
 	for ( std::size_t t = 0; t <= 10000; ++t ) {
 		computeLbmStep();
 
-		if ( t % 100 == 0 ) {
+		if ( t % 1000 == 0 ) {
 			std::cout << ".";
 			std::cout.flush();
 			fluid.writeAsVTK("result/data_t" + std::to_string(t) + ".vtk");
diff --git a/src/boundary_conditions.cc b/src/boundary_conditions.cc
index 3488463..08b157c 100644
--- a/src/boundary_conditions.cc
+++ b/src/boundary_conditions.cc
@@ -1,5 +1,7 @@
 #include "boundary_conditions.h"
 
+#include "lbm.h"
+
 std::pair<Vector<int>, Vector<int>> neighbors(Vector<int> v) {
 	if ( v[0] == 0 ) {
 		return {
@@ -27,6 +29,21 @@ void computeWallCell(DataCellBuffer& pop, Vector<std::size_t> cell, Vector<int>
 	pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB);
 }
 
+void computeMovingWallCell(DataCellBuffer& pop, Vector<std::size_t> cell, Vector<int> normal, Vector<double> u) {
+	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) - (6. * weight.get(-neighborA) * rho * (-neighborA * u));
+	pop.curr(cell).get(normal   ) = pop.curr(cell).get(-normal   ) - (6. * weight.get(-normal)    * rho * (-normal    * u));
+	pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB) - (6. * weight.get(-neighborB) * rho * (-neighborB * u));
+}
+
 void computeZouHeVelocityWallCell(DataCellBuffer& pop, Vector<std::size_t> cell, Vector<int> normal, double vX) {
 	const auto [neighborA, neighborB] = neighbors(normal);
 
diff --git a/src/boundary_conditions.h b/src/boundary_conditions.h
index 6b7de19..5f528e1 100644
--- a/src/boundary_conditions.h
+++ b/src/boundary_conditions.h
@@ -5,4 +5,6 @@
 
 void computeWallCell(DataCellBuffer& pop, Vector<std::size_t> cell, Vector<int> normal);
 
+void computeMovingWallCell(DataCellBuffer& pop, Vector<std::size_t> cell, Vector<int> normal, Vector<double> u);
+
 void computeZouHeVelocityWallCell(DataCellBuffer& pop, Vector<std::size_t> cell, Vector<int> normal, double vX);
diff --git a/src/vector.h b/src/vector.h
index 5e105c4..839707f 100644
--- a/src/vector.h
+++ b/src/vector.h
@@ -53,6 +53,11 @@ Vector<T> operator*(T scalar, const Vector<T>& v) {
 }
 
 template <typename T, typename W>
+decltype(T{}*W{}) operator*(const Vector<T>& a, const Vector<W>& b) {
+	return a[0]*b[0] + a[1]*b[1];
+}
+
+template <typename T, typename W>
 Vector<T> operator-(const Vector<T>& a, const Vector<W>& b) {
 	return Vector<T>{
 		a[0] - b[0],