Lid driven cavity with bounce-back, Zou He boundary conditions

1 files changed, 61 insertions, 72 deletions

diff --git a/boltz.cc b/boltz.cc
index 29f8518..3796b87 100644
--- a/boltz.cc
+++ b/boltz.cc
@@ -27,7 +27,6 @@ struct DataCell {
 };
 
 using Velocity = Vector<double>;
-using Force = Vector<double>;
 using Density = double;
 
 constexpr DataCell weight{
@@ -93,13 +92,12 @@ class CellBuffer {
 		}
 };
 
-constexpr std::size_t dimX = 128;
-constexpr std::size_t dimY = 128;
+constexpr std::size_t dimX = 256;
+constexpr std::size_t dimY = 256;
 
-constexpr double lidVelocityX = 0.05;
-constexpr double reynolds     = 1000;
-constexpr double nu           = lidVelocityX * dimX / reynolds;
-constexpr double tau          = 3*nu + 0.5;
+constexpr double tau          = 2./3.;
+constexpr int    reynolds     = 1000;
+constexpr double lidVelocityX = reynolds * (2.*tau - 1.) / (6.*(dimX - 1.));
 
 constexpr double omega = 1. / tau;
 
@@ -107,7 +105,6 @@ CellBuffer pop(dimX, dimY);
 
 Density  density [dimX][dimY];
 Velocity velocity[dimX][dimY];
-Force    force   [dimX][dimY];
 
 void computeFluidCell(std::size_t x, std::size_t y) {
 	// compute equilibrium
@@ -115,9 +112,21 @@ void computeFluidCell(std::size_t x, std::size_t y) {
 	eq.equilibrize(density[x][y], velocity[x][y]);
 
 	// collide (BGK, relax towards equilibrium) & stream
-	for ( int i = -1; i <= 1; ++i ) {
-		for ( int j = -1; j <= 1; ++j ) {
-			pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j) + omega * (eq.get(i,j) - pop.prev(x,y).get(i,j));
+	if ( x != 0 && x != dimX - 1 && y != 0 && y != dimY -1 ) {
+		for ( int i = -1; i <= 1; ++i ) {
+			for ( int j = -1; j <= 1; ++j ) {
+				pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j) + omega * (eq.get(i,j) - pop.prev(x,y).get(i,j));
+			}
+		}
+	} else {
+		// partial collide and stream for boundary cells,
+		// remaining populations to be determined by boundary conditions
+		for ( int i = -1; i <= 1; ++i ) {
+			for ( int j = -1; j <= 1; ++j ) {
+				if ( x+i >= 0 && x+i <= dimX - 1 && y+j >= 0 && y+j <= dimY - 1 ) {
+					pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j) + omega * (eq.get(i,j) - pop.prev(x,y).get(i,j));
+				}
+			}
 		}
 	}
 }
@@ -144,78 +153,59 @@ std::pair<Vector<int>, Vector<int>> neighbors(Vector<int> v) {
 void computeWallCell(Vector<std::size_t> cell, Vector<int> normal) {
 	const auto [neighborA, neighborB] = neighbors(normal);
 
-	pop.curr(cell).get(neighborA) = pop.curr(cell-neighborA).get(-neighborA);
-	pop.curr(cell).get(normal   ) = pop.curr(cell-normal   ).get(-normal   );
-	pop.curr(cell).get(neighborB) = pop.curr(cell-neighborB).get(-neighborB);
+	pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA);
+	pop.curr(cell).get(normal   ) = pop.curr(cell).get(-normal   );
+	pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB);
 }
 
 void computeSimpleVelocityWallCell(Vector<std::size_t> cell, Vector<int> normal, double vX) {
 	const auto [neighborA, neighborB] = neighbors(normal);
 
-	pop.curr(cell).get(neighborA) = pop.curr(cell-neighborA).get(-neighborA) - 6 * weight.get( 1, 1) * density[cell[0]][cell[1]] * vX;
-	pop.curr(cell).get(normal   ) = pop.curr(cell-normal   ).get(-normal   );
-	pop.curr(cell).get(neighborB) = pop.curr(cell-neighborB).get(-neighborB) + 6 * weight.get(-1, 1) * density[cell[0]][cell[0]] * vX;
+	pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA) - 6 * weight.get( 1, 1) * density[cell[0]][cell[1]] * vX;
+	pop.curr(cell).get(normal   ) = pop.curr(cell).get(-normal   );
+	pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB) + 6 * weight.get(-1, 1) * density[cell[0]][cell[0]] * vX;
 }
 
 void computeZouHeVelocityWallCell(Vector<std::size_t> cell, Vector<int> normal, double vX) {
 	const auto [neighborA, neighborB] = neighbors(normal);
 
-	const double rho = pop.curr(cell-neighborA).get(-1,0) + pop.curr(cell-normal).get(0,0) + pop.curr(cell-neighborB).get(1,0)
+	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-neighborA).get(-neighborA) +
-			pop.curr(cell-normal   ).get(-normal   ) +
-			pop.curr(cell-neighborB).get(-neighborB)
+			pop.curr(cell).get(-neighborA) +
+			pop.curr(cell).get(-normal   ) +
+			pop.curr(cell).get(-neighborB)
 		);
-	const double hor = pop.curr(cell-neighborA).get(-1,0) - pop.curr(cell-neighborB).get(1,0);
 
-	pop.curr(cell).get(neighborA) = pop.curr(cell-neighborA).get(-neighborA) + 0.5*hor - 0.5*vX*rho;
-	pop.curr(cell).get(normal   ) = pop.curr(cell-normal   ).get(-normal   );
-	pop.curr(cell).get(neighborB) = pop.curr(cell-neighborB).get(-neighborB) - 0.5*hor + 0.5*vX*rho;
+	pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA) + 0.5*(pop.curr(cell).get( 1,0) - pop.curr(cell).get(-1,0) - vX*rho);
+	pop.curr(cell).get(normal   ) = pop.curr(cell).get(-normal   );
+	pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB) + 0.5*(pop.curr(cell).get(-1,0) - pop.curr(cell).get( 1,0) + vX*rho);
 }
 
 void computeLbmStep(std::size_t t) {
 	pop.swap();
 
-	for ( std::size_t x = 1; x < dimX - 1; ++x ) {
-		for ( std::size_t y = 1; y < dimY - 1; ++y ) {
-			//if ( x <= 20 || x >= 40 || y <= 20 || y >= 80 ) {
-				computeFluidCell(x, y);
-			//}
+	for ( std::size_t x = 1; x < dimX-1; ++x ) {
+		for ( std::size_t y = 0; y < dimY; ++y ) {
+			computeFluidCell(x, y);
 		}
 	}
 
-	// obstacle
-	/*{
-		for ( std::size_t x = 21; x < 40; ++x ) {
-			computeWallCell(x, 80, { 0, 1});
-			computeWallCell(x, 20, { 0,-1});
-		}
-		for ( std::size_t y = 21; y < 80; ++y ) {
-			computeWallCell(40, y, { 1, 0});
-			computeWallCell(20, y, {-1, 0});
-		}
-
-		computeWallCell(40, 80, { 1, 1});
-		computeWallCell(40, 20, { 1,-1});
-		computeWallCell(20, 80, {-1, 1});
-		computeWallCell(20, 20, {-1,-1});
-	}*/
-
 	// straight wall cell bounce back
-	for ( std::size_t x = 2; x < dimX - 2; ++x ) {
-		computeWallCell({x, 1}, { 0, 1});
-		computeZouHeVelocityWallCell({x, dimY-2}, { 0,-1}, lidVelocityX);
+	for ( std::size_t x = 0; x < dimX; ++x ) {
+		computeWallCell({x, 0}, { 0, 1});
+		computeZouHeVelocityWallCell({x, dimY-1}, { 0,-1}, lidVelocityX);
+		//computeSimpleVelocityWallCell({x, dimY-1}, { 0,-1}, lidVelocityX);
 	}
-	for ( std::size_t y = 2; y < dimY - 2; ++y ) {
-		computeWallCell({1,      y}, { 1, 0});
-		computeWallCell({dimX-2, y}, {-1, 0});
+	for ( std::size_t y = 1; y < dimY-1; ++y ) {
+		computeWallCell({0,      y}, { 1, 0});
+		computeWallCell({dimX-1, y}, {-1, 0});
 	}
 
 	// edge wall cell bounce back
-	computeWallCell({1,      1     }, { 1, 1});
-	computeWallCell({1,      dimY-2}, { 1,-1});
-	computeWallCell({dimX-2, 1     }, {-1, 1});
-	computeWallCell({dimX-2, dimY-2}, {-1,-1});
+	//computeWallCell({0,      dimY-1}, { 1,-1});
+	//computeWallCell({dimX-1, dimY-1}, {-1,-1});
+	//computeWallCell({0,      0     }, { 1, 1});
+	//computeWallCell({dimX-1, 0     }, {-1, 1});
 
 	// update density, velocity field
 	for ( std::size_t x = 0; x < dimX; ++x ) {
@@ -229,9 +219,8 @@ void computeLbmStep(std::size_t t) {
 void init() {
 	for ( std::size_t x = 0; x < dimX; ++x ) {
 		for ( std::size_t y = 0; y < dimY; ++y ) {
-			density[x][y] = 1.0;
 			velocity[x][y] = { 0.0, 0.0 };
-			force[x][y] = { 0.0, 0.0 };
+			density[x][y] = 1.0;
 
 			pop.curr(x,y).equilibrize(density[x][y], velocity[x][y]);
 			pop.prev(x,y).equilibrize(density[x][y], velocity[x][y]);
@@ -247,33 +236,33 @@ void writeCurrentStateAsVTK(int time) {
 	fout << "lbm_output\n";
 	fout << "ASCII\n";
 	fout << "DATASET RECTILINEAR_GRID\n";
-	fout << "DIMENSIONS " << dimX - 2 << " " << dimY - 2 << " 1" << "\n";
+	fout << "DIMENSIONS " << dimX << " " << dimY << " 1" << "\n";
 
-	fout << "X_COORDINATES " << dimX - 2 << " float\n";
-	for( std::size_t x = 1; x < dimX - 1; ++x ) {
+	fout << "X_COORDINATES " << dimX << " float\n";
+	for( std::size_t x = 0; x < dimX; ++x ) {
 		fout << x << " ";
 	}
 
-	fout << "\nY_COORDINATES " << dimY - 2 << " float\n";
-	for( std::size_t y = 1; y < dimY - 1; ++y ) {
+	fout << "\nY_COORDINATES " << dimY << " float\n";
+	for( std::size_t y = 0; y < dimY; ++y ) {
 		fout << y << " ";
 	}
 
 	fout << "\nZ_COORDINATES " << 1 << " float\n";
 	fout << 0 << "\n";
-	fout << "POINT_DATA " << (dimX - 2) * (dimY - 2) << "\n";
+	fout << "POINT_DATA " << dimX * dimY << "\n";
 
 	fout << "VECTORS velocity float\n";
-	for ( std::size_t y = 1; y < dimY - 1; ++y ) {
-		for ( std::size_t x = 1; x < dimX - 1; ++x ) {
+	for ( std::size_t y = 0; y < dimY; ++y ) {
+		for ( std::size_t x = 0; x < dimX; ++x ) {
 			fout << velocity[x][y][0] << " " << velocity[x][y][1] << " 0\n";
 		}
 	}
 
 	fout << "SCALARS density float 1\n";
 	fout << "LOOKUP_TABLE default\n";
-	for ( std::size_t y = 1; y < dimY - 1; ++y ) {
-		for ( std::size_t x = 1; x < dimX - 1; ++x ) {
+	for ( std::size_t y = 0; y < dimY; ++y ) {
+		for ( std::size_t x = 0; x < dimX; ++x ) {
 			fout << density[x][y] << "\n";
 		}
 	}
@@ -284,13 +273,13 @@ void writeCurrentStateAsVTK(int time) {
 int main() {
 	init();
 
-	std::cout << "u: " << lidVelocityX << std::endl;
 	std::cout << "tau: " << tau << std::endl;
+	std::cout << "u: " << lidVelocityX << std::endl;
 
-	for ( std::size_t t = 0; t < 20000; ++t ) {
+	for ( std::size_t t = 0; t <= 40000; ++t ) {
 		computeLbmStep(t);
 
-		if ( t % 1000 == 0 ) {
+		if ( t < 100 || t % 1000 == 0 ) {
 			std::cout << "Writing " << t << "." << std::endl;
 			writeCurrentStateAsVTK(t);
 		}