Extract public version of SweepLB

1 files changed, 87 insertions, 0 deletions

diff --git a/channel.cpp b/channel.cpp
new file mode 100644
index 0000000..fddb9a1
--- /dev/null
+++ b/channel.cpp
@@ -0,0 +1,87 @@
+#include "lattice.h"
+
+#include "LLBM/collide.h"
+#include "LLBM/initialize.h"
+#include "LLBM/bounce_back.h"
+#include "LLBM/bounce_back_moving_wall.h"
+#include "LLBM/equilibrium_density_wall.h"
+#include "LLBM/equilibrium_velocity_wall.h"
+
+#include <chrono>
+#include <iostream>
+
+#include "pattern/all.h"
+
+using T = SWEEPLB_PRECISION;
+using PATTERN = pattern::SWEEPLB_PATTERN<T>;
+
+void simulate(Cuboid cuboid, std::size_t nStep) {
+	const int nThread = omp_get_max_threads();
+
+	Lattice<PATTERN> lattice(cuboid);
+
+	LatticeMask<T> bulk_mask(cuboid.volume());
+	LatticeMask<T> wall_mask(cuboid.volume());
+	LatticeMask<T> inflow_mask(cuboid.volume());
+	LatticeMask<T> outflow_mask(cuboid.volume());
+
+	cuboid.traverse([&](int iX, int iY, int iZ, std::size_t iCell) {
+		if (   iY == 0 || iY == cuboid[1]-1
+		    || iZ == 0 || iZ == cuboid[2]-1) {
+			wall_mask.set(iCell, true);
+		} else if (iX == 0) {
+			inflow_mask.set(iCell, true);
+		} else if (iX == cuboid[0]-1) {
+			outflow_mask.set(iCell, true);
+		} else {
+			bulk_mask.set(iCell, true);
+		}
+	});
+
+	bulk_mask.serialize();
+	wall_mask.serialize();
+	inflow_mask.serialize();
+	outflow_mask.serialize();
+
+	T tau = 0.56;
+
+	T u_inflow  = 0.05;
+	T d_outflow = 1.;
+
+	for (std::size_t iStep = 0; iStep < 100; ++iStep) {
+		lattice.apply(Operator(BgkCollideO(), bulk_mask, tau),
+		              Operator(BounceBackO(), wall_mask),
+		              Operator(EquilibriumVelocityWallO(), inflow_mask, u_inflow, WallNormal<1,0,0>()),
+		              Operator(EquilibriumDensityWallO(), outflow_mask, d_outflow, WallNormal<-1,0,0>()));
+		lattice.stream();
+	}
+
+	auto start = std::chrono::steady_clock::now();
+
+	for (std::size_t iStep = 0; iStep < nStep; ++iStep) {
+		lattice.apply(Operator(BgkCollideO(), bulk_mask, tau),
+		              Operator(BounceBackO(), wall_mask),
+		              Operator(EquilibriumVelocityWallO(), inflow_mask, u_inflow, WallNormal<1,0,0>()),
+		              Operator(EquilibriumDensityWallO(), outflow_mask, d_outflow, WallNormal<-1,0,0>()));
+		lattice.stream();
+	}
+
+	auto duration = std::chrono::duration_cast<std::chrono::duration<double>>(
+	  std::chrono::steady_clock::now() - start);
+
+	std::cout << cuboid[0] << ", " << cuboid[1] << ", " << cuboid[2]
+	          << ", " << nStep
+	          << ", " << nThread
+	          << ", " << (nStep * lattice.volume()) / (1e6 * duration.count())
+	          << std::endl;
+
+	lattice.write_momenta(bulk_mask, "result.vtk");
+}
+
+int main(int argc, char* argv[]) {
+	const std::size_t nX    = atoi(argv[1]);
+	const std::size_t nY    = atoi(argv[2]);
+	const std::size_t nZ    = atoi(argv[3]);
+	const std::size_t steps = atoi(argv[4]);
+	simulate({ nX, nY, nZ }, steps);
+}