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.

1 files changed, 17 insertions, 0 deletions

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);