Extract public version of SweepLB

9 files changed, 648 insertions, 0 deletions

diff --git a/src/LLBM/bounce_back.h b/src/LLBM/bounce_back.h
new file mode 100644
index 0000000..b710ee4
--- /dev/null
+++ b/src/LLBM/bounce_back.h
@@ -0,0 +1,30 @@
+#pragma once
+
+#include "concepts.h"
+
+struct BounceBackO {
+
+template <concepts::Arithmetic V>
+static void apply(V f_curr[19], V f_next[19]) {
+	f_next[0] = f_curr[18];
+	f_next[1] = f_curr[17];
+	f_next[2] = f_curr[16];
+	f_next[3] = f_curr[15];
+	f_next[4] = f_curr[14];
+	f_next[5] = f_curr[13];
+	f_next[6] = f_curr[12];
+	f_next[7] = f_curr[11];
+	f_next[8] = f_curr[10];
+	f_next[9] = f_curr[9];
+	f_next[10] = f_curr[8];
+	f_next[11] = f_curr[7];
+	f_next[12] = f_curr[6];
+	f_next[13] = f_curr[5];
+	f_next[14] = f_curr[4];
+	f_next[15] = f_curr[3];
+	f_next[16] = f_curr[2];
+	f_next[17] = f_curr[1];
+	f_next[18] = f_curr[0];
+}
+
+};
diff --git a/src/LLBM/bounce_back_moving_wall.h b/src/LLBM/bounce_back_moving_wall.h
new file mode 100644
index 0000000..c00a17e
--- /dev/null
+++ b/src/LLBM/bounce_back_moving_wall.h
@@ -0,0 +1,30 @@
+#pragma once
+
+#include "concepts.h"
+
+struct BounceBackMovingWallO {
+
+template <concepts::Arithmetic V, concepts::Arithmetic T>
+static void apply(V f_curr[19], V f_next[19], T u[3]) {
+	f_next[0] = f_curr[18] + T{0.166666666666667}*u[1] + T{0.166666666666667}*u[2];
+	f_next[1] = f_curr[17] - T{0.166666666666667}*u[0] + T{0.166666666666667}*u[2];
+	f_next[2] = f_curr[16] + T{0.333333333333333}*u[2];
+	f_next[3] = f_curr[15] + T{0.166666666666667}*u[0] + T{0.166666666666667}*u[2];
+	f_next[4] = f_curr[14] - T{0.166666666666667}*u[1] + T{0.166666666666667}*u[2];
+	f_next[5] = f_curr[13] - T{0.166666666666667}*u[0] + T{0.166666666666667}*u[1];
+	f_next[6] = f_curr[12] + T{0.333333333333333}*u[1];
+	f_next[7] = f_curr[11] + T{0.166666666666667}*u[0] + T{0.166666666666667}*u[1];
+	f_next[8] = f_curr[10] - T{0.333333333333333}*u[0];
+	f_next[9] = f_curr[9];
+	f_next[10] = f_curr[8] + T{0.333333333333333}*u[0];
+	f_next[11] = f_curr[7] - T{0.166666666666667}*u[0] - T{0.166666666666667}*u[1];
+	f_next[12] = f_curr[6] - T{0.333333333333333}*u[1];
+	f_next[13] = f_curr[5] + T{0.166666666666667}*u[0] - T{0.166666666666667}*u[1];
+	f_next[14] = f_curr[4] + T{0.166666666666667}*u[1] - T{0.166666666666667}*u[2];
+	f_next[15] = f_curr[3] - T{0.166666666666667}*u[0] - T{0.166666666666667}*u[2];
+	f_next[16] = f_curr[2] - T{0.333333333333333}*u[2];
+	f_next[17] = f_curr[1] + T{0.166666666666667}*u[0] - T{0.166666666666667}*u[2];
+	f_next[18] = f_curr[0] - T{0.166666666666667}*u[1] - T{0.166666666666667}*u[2];
+}
+
+};
diff --git a/src/LLBM/collect_moments.h b/src/LLBM/collect_moments.h
new file mode 100644
index 0000000..afd5d59
--- /dev/null
+++ b/src/LLBM/collect_moments.h
@@ -0,0 +1,25 @@
+#pragma once
+
+#include "concepts.h"
+
+struct CollectMomentsF {
+
+static constexpr unsigned size = 4;
+
+template <concepts::Arithmetic V>
+static void observe(V f_curr[19], V* momenta) {
+	V m0 = f_curr[10] + f_curr[13] + f_curr[17];
+	V m1 = f_curr[14] + f_curr[5] + f_curr[6];
+	V m2 = f_curr[1] + f_curr[2] + f_curr[4];
+	V m3 = m0 + m1 + m2 + f_curr[0] + f_curr[11] + f_curr[12] + f_curr[15] + f_curr[16] + f_curr[18] + f_curr[3] + f_curr[7] + f_curr[8] + f_curr[9];
+	V m4 = -f_curr[11] + f_curr[7];
+	V m5 = -f_curr[15] + f_curr[3];
+	V m6 = V{1}/m3;
+	V m7 = f_curr[0] - f_curr[18];
+	momenta[0] = m3;
+	momenta[1] = m6*(m0 + m4 + m5 - f_curr[1] - f_curr[5] - f_curr[8]);
+	momenta[2] = m6*(m1 + m4 + m7 - f_curr[12] - f_curr[13] - f_curr[4]);
+	momenta[3] = m6*(m2 + m5 + m7 - f_curr[14] - f_curr[16] - f_curr[17]);
+}
+
+};
diff --git a/src/LLBM/collide.h b/src/LLBM/collide.h
new file mode 100644
index 0000000..c08aade
--- /dev/null
+++ b/src/LLBM/collide.h
@@ -0,0 +1,89 @@
+#pragma once
+
+#include "concepts.h"
+
+struct BgkCollideO {
+
+template <concepts::Arithmetic V, concepts::Arithmetic T>
+static void apply(V f_curr[19], V f_next[19], T tau) {
+	V rho;
+	V u[3];
+	V m0 = f_curr[10] + f_curr[13] + f_curr[17];
+	V m1 = f_curr[14] + f_curr[5] + f_curr[6];
+	V m2 = f_curr[1] + f_curr[2] + f_curr[4];
+	V m3 = m0 + m1 + m2 + f_curr[0] + f_curr[11] + f_curr[12] + f_curr[15] + f_curr[16] + f_curr[18] + f_curr[3] + f_curr[7] + f_curr[8] + f_curr[9];
+	V m4 = -f_curr[11] + f_curr[7];
+	V m5 = -f_curr[15] + f_curr[3];
+	V m6 = T{1}/m3;
+	V m7 = f_curr[0] - f_curr[18];
+	rho = m3;
+	u[0] = m6*(m0 + m4 + m5 - f_curr[1] - f_curr[5] - f_curr[8]);
+	u[1] = m6*(m1 + m4 + m7 - f_curr[12] - f_curr[13] - f_curr[4]);
+	u[2] = m6*(m2 + m5 + m7 - f_curr[14] - f_curr[16] - f_curr[17]);
+	V x0 = T{1}/tau;
+	V x1 = T{0.0138888888888889}*x0;
+	V x2 = u[1] + u[2];
+	V x3 = T{9.00000000000000}*(x2*x2);
+	V x4 = T{6.00000000000000}*u[2];
+	V x5 = u[1]*u[1];
+	V x6 = T{3.00000000000000}*x5;
+	V x7 = -x6;
+	V x8 = u[0]*u[0];
+	V x9 = T{3.00000000000000}*x8;
+	V x10 = T{2.00000000000000} - x9;
+	V x11 = x10 + x7;
+	V x12 = x11 + x4;
+	V x13 = u[2]*u[2];
+	V x14 = T{3.00000000000000}*x13;
+	V x15 = -x14;
+	V x16 = T{6.00000000000000}*u[1];
+	V x17 = x15 + x16;
+	V x18 = T{6.00000000000000}*u[0];
+	V x19 = -u[0];
+	V x20 = x19 + u[2];
+	V x21 = x14 + x6 + T{-2.00000000000000};
+	V x22 = x21 + x9;
+	V x23 = x22 - x4;
+	V x24 = T{0.0277777777777778}*x0;
+	V x25 = T{6.00000000000000}*x13;
+	V x26 = u[0] + u[2];
+	V x27 = T{9.00000000000000}*(x26*x26);
+	V x28 = x15 + x18;
+	V x29 = -u[1];
+	V x30 = x29 + u[2];
+	V x31 = x19 + u[1];
+	V x32 = -x16;
+	V x33 = x18 + x22;
+	V x34 = T{6.00000000000000}*x5;
+	V x35 = u[0] + u[1];
+	V x36 = T{9.00000000000000}*(x35*x35);
+	V x37 = T{6.00000000000000}*x8;
+	V x38 = x9 + T{-2.00000000000000};
+	V x39 = x29 + u[0];
+	V x40 = -x18 + x22;
+	V x41 = -u[2];
+	V x42 = x41 + u[1];
+	V x43 = x22 + x4;
+	V x44 = x41 + u[0];
+	f_next[0] = x1*(rho*(x12 + x17 + x3) - T{72.0000000000000}*f_curr[0]) + f_curr[0];
+	f_next[1] = -x1*(rho*(x18 + x23 - T{9.00000000000000}*x20*x20) + T{72.0000000000000}*f_curr[1]) + f_curr[1];
+	f_next[2] = x24*(rho*(x12 + x25) - T{36.0000000000000}*f_curr[2]) + f_curr[2];
+	f_next[3] = x1*(rho*(x12 + x27 + x28) - T{72.0000000000000}*f_curr[3]) + f_curr[3];
+	f_next[4] = -x1*(rho*(x16 + x23 - T{9.00000000000000}*x30*x30) + T{72.0000000000000}*f_curr[4]) + f_curr[4];
+	f_next[5] = -x1*(rho*(x32 + x33 - T{9.00000000000000}*x31*x31) + T{72.0000000000000}*f_curr[5]) + f_curr[5];
+	f_next[6] = x24*(rho*(x10 + x17 + x34) - T{36.0000000000000}*f_curr[6]) + f_curr[6];
+	f_next[7] = x1*(rho*(x11 + x17 + x18 + x36) - T{72.0000000000000}*f_curr[7]) + f_curr[7];
+	f_next[8] = -x24*(rho*(x18 + x21 - x37) + T{36.0000000000000}*f_curr[8]) + f_curr[8];
+	f_next[9] = -T{0.166666666666667}*x0*(rho*x22 + T{6.00000000000000}*f_curr[9]) + f_curr[9];
+	f_next[10] = x24*(rho*(x28 + x37 + x7 + T{2.00000000000000}) - T{36.0000000000000}*f_curr[10]) + f_curr[10];
+	f_next[11] = -x1*(rho*(x16 + x33 - x36) + T{72.0000000000000}*f_curr[11]) + f_curr[11];
+	f_next[12] = -x24*(rho*(x14 + x16 - x34 + x38) + T{36.0000000000000}*f_curr[12]) + f_curr[12];
+	f_next[13] = -x1*(rho*(x16 + x40 - T{9.00000000000000}*x39*x39) + T{72.0000000000000}*f_curr[13]) + f_curr[13];
+	f_next[14] = -x1*(rho*(x32 + x43 - T{9.00000000000000}*x42*x42) + T{72.0000000000000}*f_curr[14]) + f_curr[14];
+	f_next[15] = -x1*(rho*(-x27 + x33 + x4) + T{72.0000000000000}*f_curr[15]) + f_curr[15];
+	f_next[16] = -x24*(rho*(-x25 + x38 + x4 + x6) + T{36.0000000000000}*f_curr[16]) + f_curr[16];
+	f_next[17] = -x1*(rho*(x4 + x40 - T{9.00000000000000}*x44*x44) + T{72.0000000000000}*f_curr[17]) + f_curr[17];
+	f_next[18] = -x1*(rho*(x16 - x3 + x43) + T{72.0000000000000}*f_curr[18]) + f_curr[18];
+}
+
+};
diff --git a/src/LLBM/equilibrium_density_wall.h b/src/LLBM/equilibrium_density_wall.h
new file mode 100644
index 0000000..598d59d
--- /dev/null
+++ b/src/LLBM/equilibrium_density_wall.h
@@ -0,0 +1,144 @@
+#pragma once
+
+#include "wall.h"
+#include "concepts.h"
+
+struct EquilibriumDensityWallO {
+
+template <concepts::Arithmetic V, concepts::Arithmetic T>
+static void apply(V f_curr[19], V f_next[19], T rho_w, WallNormal<1,0,0>) {
+	V u_w = (rho_w - f_curr[0] - T{2.00000000000000}*f_curr[11] - f_curr[12] - f_curr[14] - T{2.00000000000000}*f_curr[15] - f_curr[16] - f_curr[18] - T{2.00000000000000}*f_curr[1] - f_curr[2] - f_curr[4] - T{2.00000000000000}*f_curr[5] - f_curr[6] - T{2.00000000000000}*f_curr[8] - f_curr[9])/rho_w;
+	V rho = rho_w;
+	V u[3] { u_w, 0., 0. };
+	V e0 = T{0.0277777777777778}*rho;
+	V e1 = u[1] + u[2];
+	V e2 = T{4.50000000000000}*(e1*e1);
+	V e3 = T{3.00000000000000}*u[2];
+	V e4 = u[1]*u[1];
+	V e5 = T{1.50000000000000}*e4;
+	V e6 = -e5;
+	V e7 = u[0]*u[0];
+	V e8 = T{1.50000000000000}*e7;
+	V e9 = T{1.00000000000000} - e8;
+	V e10 = e6 + e9;
+	V e11 = e10 + e3;
+	V e12 = T{3.00000000000000}*u[1];
+	V e13 = u[2]*u[2];
+	V e14 = T{1.50000000000000}*e13;
+	V e15 = -e14;
+	V e16 = e12 + e15;
+	V e17 = T{3.00000000000000}*u[0];
+	V e18 = -u[2];
+	V e19 = e18 + u[0];
+	V e20 = -T{4.50000000000000}*e19*e19;
+	V e21 = e14 + e5 + T{-1.00000000000000};
+	V e22 = e21 + e8;
+	V e23 = e22 - e3;
+	V e24 = T{0.0555555555555556}*rho;
+	V e25 = T{3.00000000000000}*e13;
+	V e26 = u[0] + u[2];
+	V e27 = T{4.50000000000000}*(e26*e26);
+	V e28 = e15 + e17;
+	V e29 = e18 + u[1];
+	V e30 = -T{4.50000000000000}*e29*e29;
+	V e31 = -e12;
+	V e32 = u[0] - u[1];
+	V e33 = -T{4.50000000000000}*e32*e32;
+	V e34 = e17 + e22;
+	V e35 = T{3.00000000000000}*e4;
+	V e36 = u[0] + u[1];
+	V e37 = T{4.50000000000000}*(e36*e36);
+	V e38 = T{3.00000000000000}*e7;
+	V e39 = e8 + T{-1.00000000000000};
+	V e40 = -e17 + e22;
+	V e41 = e22 + e3;
+	f_next[0] = e0*(e11 + e16 + e2);
+	f_next[1] = -e0*(e17 + e20 + e23);
+	f_next[2] = e24*(e11 + e25);
+	f_next[3] = e0*(e11 + e27 + e28);
+	f_next[4] = -e0*(e12 + e23 + e30);
+	f_next[5] = -e0*(e31 + e33 + e34);
+	f_next[6] = e24*(e16 + e35 + e9);
+	f_next[7] = e0*(e10 + e16 + e17 + e37);
+	f_next[8] = -e24*(e17 + e21 - e38);
+	f_next[9] = -T{0.333333333333333}*e22*rho;
+	f_next[10] = e24*(e28 + e38 + e6 + T{1.00000000000000});
+	f_next[11] = -e0*(e12 + e34 - e37);
+	f_next[12] = -e24*(e12 + e14 - e35 + e39);
+	f_next[13] = -e0*(e12 + e33 + e40);
+	f_next[14] = -e0*(e30 + e31 + e41);
+	f_next[15] = -e0*(-e27 + e3 + e34);
+	f_next[16] = -e24*(-e25 + e3 + e39 + e5);
+	f_next[17] = -e0*(e20 + e3 + e40);
+	f_next[18] = -e0*(e12 - e2 + e41);
+}
+
+template <concepts::Arithmetic V, concepts::Arithmetic T>
+static void apply(V f_curr[19], V f_next[19], T rho_w, WallNormal<-1,0,0>) {
+	V u_w = (-rho_w + f_curr[0] + T{2.00000000000000}*f_curr[10] + f_curr[12] + T{2.00000000000000}*f_curr[13] + f_curr[14] + f_curr[16] + T{2.00000000000000}*f_curr[17] + f_curr[18] + f_curr[2] + T{2.00000000000000}*f_curr[3] + f_curr[4] + f_curr[6] + T{2.00000000000000}*f_curr[7] + f_curr[9])/rho_w;
+	V rho = rho_w;
+	V u[3] { u_w, 0., 0. };
+	V e0 = T{0.0277777777777778}*rho;
+	V e1 = u[1] + u[2];
+	V e2 = T{4.50000000000000}*(e1*e1);
+	V e3 = T{3.00000000000000}*u[2];
+	V e4 = u[1]*u[1];
+	V e5 = T{1.50000000000000}*e4;
+	V e6 = -e5;
+	V e7 = u[0]*u[0];
+	V e8 = T{1.50000000000000}*e7;
+	V e9 = T{1.00000000000000} - e8;
+	V e10 = e6 + e9;
+	V e11 = e10 + e3;
+	V e12 = T{3.00000000000000}*u[1];
+	V e13 = u[2]*u[2];
+	V e14 = T{1.50000000000000}*e13;
+	V e15 = -e14;
+	V e16 = e12 + e15;
+	V e17 = T{3.00000000000000}*u[0];
+	V e18 = -u[2];
+	V e19 = e18 + u[0];
+	V e20 = -T{4.50000000000000}*e19*e19;
+	V e21 = e14 + e5 + T{-1.00000000000000};
+	V e22 = e21 + e8;
+	V e23 = e22 - e3;
+	V e24 = T{0.0555555555555556}*rho;
+	V e25 = T{3.00000000000000}*e13;
+	V e26 = u[0] + u[2];
+	V e27 = T{4.50000000000000}*(e26*e26);
+	V e28 = e15 + e17;
+	V e29 = e18 + u[1];
+	V e30 = -T{4.50000000000000}*e29*e29;
+	V e31 = -e12;
+	V e32 = u[0] - u[1];
+	V e33 = -T{4.50000000000000}*e32*e32;
+	V e34 = e17 + e22;
+	V e35 = T{3.00000000000000}*e4;
+	V e36 = u[0] + u[1];
+	V e37 = T{4.50000000000000}*(e36*e36);
+	V e38 = T{3.00000000000000}*e7;
+	V e39 = e8 + T{-1.00000000000000};
+	V e40 = -e17 + e22;
+	V e41 = e22 + e3;
+	f_next[0] = e0*(e11 + e16 + e2);
+	f_next[1] = -e0*(e17 + e20 + e23);
+	f_next[2] = e24*(e11 + e25);
+	f_next[3] = e0*(e11 + e27 + e28);
+	f_next[4] = -e0*(e12 + e23 + e30);
+	f_next[5] = -e0*(e31 + e33 + e34);
+	f_next[6] = e24*(e16 + e35 + e9);
+	f_next[7] = e0*(e10 + e16 + e17 + e37);
+	f_next[8] = -e24*(e17 + e21 - e38);
+	f_next[9] = -T{0.333333333333333}*e22*rho;
+	f_next[10] = e24*(e28 + e38 + e6 + T{1.00000000000000});
+	f_next[11] = -e0*(e12 + e34 - e37);
+	f_next[12] = -e24*(e12 + e14 - e35 + e39);
+	f_next[13] = -e0*(e12 + e33 + e40);
+	f_next[14] = -e0*(e30 + e31 + e41);
+	f_next[15] = -e0*(-e27 + e3 + e34);
+	f_next[16] = -e24*(-e25 + e3 + e39 + e5);
+	f_next[17] = -e0*(e20 + e3 + e40);
+	f_next[18] = -e0*(e12 - e2 + e41);
+}
+
+};
diff --git a/src/LLBM/equilibrium_velocity_wall.h b/src/LLBM/equilibrium_velocity_wall.h
new file mode 100644
index 0000000..2d68e34
--- /dev/null
+++ b/src/LLBM/equilibrium_velocity_wall.h
@@ -0,0 +1,144 @@
+#pragma once
+
+#include "wall.h"
+#include "concepts.h"
+
+struct EquilibriumVelocityWallO {
+
+template <concepts::Arithmetic V, concepts::Arithmetic T>
+static void apply(V f_curr[19], V f_next[19], T u_w, WallNormal<1,0,0>) {
+	V rho;
+	rho = -(f_curr[0] + T{2.00000000000000}*f_curr[11] + f_curr[12] + f_curr[14] + T{2.00000000000000}*f_curr[15] + f_curr[16] + f_curr[18] + T{2.00000000000000}*f_curr[1] + f_curr[2] + f_curr[4] + T{2.00000000000000}*f_curr[5] + f_curr[6] + T{2.00000000000000}*f_curr[8] + f_curr[9])/(u_w + T{-1.00000000000000});
+	V u[3] { u_w, 0., 0. };
+	V e0 = T{0.0277777777777778}*rho;
+	V e1 = u[1] + u[2];
+	V e2 = T{4.50000000000000}*(e1*e1);
+	V e3 = T{3.00000000000000}*u[2];
+	V e4 = u[1]*u[1];
+	V e5 = T{1.50000000000000}*e4;
+	V e6 = -e5;
+	V e7 = u[0]*u[0];
+	V e8 = T{1.50000000000000}*e7;
+	V e9 = T{1.00000000000000} - e8;
+	V e10 = e6 + e9;
+	V e11 = e10 + e3;
+	V e12 = T{3.00000000000000}*u[1];
+	V e13 = u[2]*u[2];
+	V e14 = T{1.50000000000000}*e13;
+	V e15 = -e14;
+	V e16 = e12 + e15;
+	V e17 = T{3.00000000000000}*u[0];
+	V e18 = -u[2];
+	V e19 = e18 + u[0];
+	V e20 = -T{4.50000000000000}*e19*e19;
+	V e21 = e14 + e5 + T{-1.00000000000000};
+	V e22 = e21 + e8;
+	V e23 = e22 - e3;
+	V e24 = T{0.0555555555555556}*rho;
+	V e25 = T{3.00000000000000}*e13;
+	V e26 = u[0] + u[2];
+	V e27 = T{4.50000000000000}*(e26*e26);
+	V e28 = e15 + e17;
+	V e29 = e18 + u[1];
+	V e30 = -T{4.50000000000000}*e29*e29;
+	V e31 = -e12;
+	V e32 = u[0] - u[1];
+	V e33 = -T{4.50000000000000}*e32*e32;
+	V e34 = e17 + e22;
+	V e35 = T{3.00000000000000}*e4;
+	V e36 = u[0] + u[1];
+	V e37 = T{4.50000000000000}*(e36*e36);
+	V e38 = T{3.00000000000000}*e7;
+	V e39 = e8 + T{-1.00000000000000};
+	V e40 = -e17 + e22;
+	V e41 = e22 + e3;
+	f_next[0] = e0*(e11 + e16 + e2);
+	f_next[1] = -e0*(e17 + e20 + e23);
+	f_next[2] = e24*(e11 + e25);
+	f_next[3] = e0*(e11 + e27 + e28);
+	f_next[4] = -e0*(e12 + e23 + e30);
+	f_next[5] = -e0*(e31 + e33 + e34);
+	f_next[6] = e24*(e16 + e35 + e9);
+	f_next[7] = e0*(e10 + e16 + e17 + e37);
+	f_next[8] = -e24*(e17 + e21 - e38);
+	f_next[9] = -T{0.333333333333333}*e22*rho;
+	f_next[10] = e24*(e28 + e38 + e6 + T{1.00000000000000});
+	f_next[11] = -e0*(e12 + e34 - e37);
+	f_next[12] = -e24*(e12 + e14 - e35 + e39);
+	f_next[13] = -e0*(e12 + e33 + e40);
+	f_next[14] = -e0*(e30 + e31 + e41);
+	f_next[15] = -e0*(-e27 + e3 + e34);
+	f_next[16] = -e24*(-e25 + e3 + e39 + e5);
+	f_next[17] = -e0*(e20 + e3 + e40);
+	f_next[18] = -e0*(e12 - e2 + e41);
+}
+
+template <concepts::Arithmetic V, concepts::Arithmetic T>
+static void apply(V f_curr[19], V f_next[19], T u_w, WallNormal<-1,0,0>) {
+	V rho;
+	rho = (f_curr[0] + T{2.00000000000000}*f_curr[10] + f_curr[12] + T{2.00000000000000}*f_curr[13] + f_curr[14] + f_curr[16] + T{2.00000000000000}*f_curr[17] + f_curr[18] + f_curr[2] + T{2.00000000000000}*f_curr[3] + f_curr[4] + f_curr[6] + T{2.00000000000000}*f_curr[7] + f_curr[9])/(u_w + T{1.00000000000000});
+	V u[3] { u_w, 0., 0. };
+	V e0 = T{0.0277777777777778}*rho;
+	V e1 = u[1] + u[2];
+	V e2 = T{4.50000000000000}*(e1*e1);
+	V e3 = T{3.00000000000000}*u[2];
+	V e4 = u[1]*u[1];
+	V e5 = T{1.50000000000000}*e4;
+	V e6 = -e5;
+	V e7 = u[0]*u[0];
+	V e8 = T{1.50000000000000}*e7;
+	V e9 = T{1.00000000000000} - e8;
+	V e10 = e6 + e9;
+	V e11 = e10 + e3;
+	V e12 = T{3.00000000000000}*u[1];
+	V e13 = u[2]*u[2];
+	V e14 = T{1.50000000000000}*e13;
+	V e15 = -e14;
+	V e16 = e12 + e15;
+	V e17 = T{3.00000000000000}*u[0];
+	V e18 = -u[2];
+	V e19 = e18 + u[0];
+	V e20 = -T{4.50000000000000}*e19*e19;
+	V e21 = e14 + e5 + T{-1.00000000000000};
+	V e22 = e21 + e8;
+	V e23 = e22 - e3;
+	V e24 = T{0.0555555555555556}*rho;
+	V e25 = T{3.00000000000000}*e13;
+	V e26 = u[0] + u[2];
+	V e27 = T{4.50000000000000}*(e26*e26);
+	V e28 = e15 + e17;
+	V e29 = e18 + u[1];
+	V e30 = -T{4.50000000000000}*e29*e29;
+	V e31 = -e12;
+	V e32 = u[0] - u[1];
+	V e33 = -T{4.50000000000000}*e32*e32;
+	V e34 = e17 + e22;
+	V e35 = T{3.00000000000000}*e4;
+	V e36 = u[0] + u[1];
+	V e37 = T{4.50000000000000}*(e36*e36);
+	V e38 = T{3.00000000000000}*e7;
+	V e39 = e8 + T{-1.00000000000000};
+	V e40 = -e17 + e22;
+	V e41 = e22 + e3;
+	f_next[0] = e0*(e11 + e16 + e2);
+	f_next[1] = -e0*(e17 + e20 + e23);
+	f_next[2] = e24*(e11 + e25);
+	f_next[3] = e0*(e11 + e27 + e28);
+	f_next[4] = -e0*(e12 + e23 + e30);
+	f_next[5] = -e0*(e31 + e33 + e34);
+	f_next[6] = e24*(e16 + e35 + e9);
+	f_next[7] = e0*(e10 + e16 + e17 + e37);
+	f_next[8] = -e24*(e17 + e21 - e38);
+	f_next[9] = -T{0.333333333333333}*e22*rho;
+	f_next[10] = e24*(e28 + e38 + e6 + T{1.00000000000000});
+	f_next[11] = -e0*(e12 + e34 - e37);
+	f_next[12] = -e24*(e12 + e14 - e35 + e39);
+	f_next[13] = -e0*(e12 + e33 + e40);
+	f_next[14] = -e0*(e30 + e31 + e41);
+	f_next[15] = -e0*(-e27 + e3 + e34);
+	f_next[16] = -e24*(-e25 + e3 + e39 + e5);
+	f_next[17] = -e0*(e20 + e3 + e40);
+	f_next[18] = -e0*(e12 - e2 + e41);
+}
+
+};
diff --git a/src/LLBM/initialize.h b/src/LLBM/initialize.h
new file mode 100644
index 0000000..b361217
--- /dev/null
+++ b/src/LLBM/initialize.h
@@ -0,0 +1,53 @@
+#pragma once
+
+#include "concepts.h"
+
+struct InitializeO {
+
+template <concepts::Arithmetic V, concepts::Arithmetic T = V>
+static void apply(V f_curr[19], V f_next[19]) {
+	f_next[0] = T{0.0277777777777778};
+	f_next[1] = T{0.0277777777777778};
+	f_next[2] = T{0.0555555555555556};
+	f_next[3] = T{0.0277777777777778};
+	f_next[4] = T{0.0277777777777778};
+	f_next[5] = T{0.0277777777777778};
+	f_next[6] = T{0.0555555555555556};
+	f_next[7] = T{0.0277777777777778};
+	f_next[8] = T{0.0555555555555556};
+	f_next[9] = T{0.333333333333333};
+	f_next[10] = T{0.0555555555555556};
+	f_next[11] = T{0.0277777777777778};
+	f_next[12] = T{0.0555555555555556};
+	f_next[13] = T{0.0277777777777778};
+	f_next[14] = T{0.0277777777777778};
+	f_next[15] = T{0.0277777777777778};
+	f_next[16] = T{0.0555555555555556};
+	f_next[17] = T{0.0277777777777778};
+	f_next[18] = T{0.0277777777777778};
+}
+
+template <concepts::Arithmetic V, concepts::Arithmetic T = V>
+static void apply(V f_curr[19], V f_next[19], T scale) {
+	f_next[0]  = T{scale * 0.0277777777777778};
+	f_next[1]  = T{scale * 0.0277777777777778};
+	f_next[2]  = T{scale * 0.0555555555555556};
+	f_next[3]  = T{scale * 0.0277777777777778};
+	f_next[4]  = T{scale * 0.0277777777777778};
+	f_next[5]  = T{scale * 0.0277777777777778};
+	f_next[6]  = T{scale * 0.0555555555555556};
+	f_next[7]  = T{scale * 0.0277777777777778};
+	f_next[8]  = T{scale * 0.0555555555555556};
+	f_next[9]  = T{scale * 0.333333333333333};
+	f_next[10] = T{scale * 0.0555555555555556};
+	f_next[11] = T{scale * 0.0277777777777778};
+	f_next[12] = T{scale * 0.0555555555555556};
+	f_next[13] = T{scale * 0.0277777777777778};
+	f_next[14] = T{scale * 0.0277777777777778};
+	f_next[15] = T{scale * 0.0277777777777778};
+	f_next[16] = T{scale * 0.0555555555555556};
+	f_next[17] = T{scale * 0.0277777777777778};
+	f_next[18] = T{scale * 0.0277777777777778};
+}
+
+};
diff --git a/src/LLBM/smagorinsky_collide.h b/src/LLBM/smagorinsky_collide.h
new file mode 100644
index 0000000..1a90598
--- /dev/null
+++ b/src/LLBM/smagorinsky_collide.h
@@ -0,0 +1,129 @@
+#pragma once
+
+#include "concepts.h"
+
+struct SmagorinskyBgkCollideO {
+
+template <concepts::Arithmetic V, concepts::Arithmetic T>
+static void apply(V f_curr[19], V f_next[19], T tau, T smagorinsky) {
+	V rho;
+	V u[3];
+	V x0 = f_curr[10] + f_curr[13] + f_curr[17];
+	V x1 = f_curr[14] + f_curr[5] + f_curr[6];
+	V x2 = f_curr[1] + f_curr[2] + f_curr[4];
+	V x3 = x0 + x1 + x2 + f_curr[0] + f_curr[11] + f_curr[12] + f_curr[15] + f_curr[16] + f_curr[18] + f_curr[3] + f_curr[7] + f_curr[8] + f_curr[9];
+	V x4 = -f_curr[11] + f_curr[7];
+	V x5 = -f_curr[15] + f_curr[3];
+	V x6 = T{1}/x3;
+	V x7 = f_curr[0] - f_curr[18];
+	V x8 = T{2.00000000000000}*tau + T{-1.00000000000000};
+	V x9 = T{72.0000000000000}*f_curr[5];
+	V x26 = T{72.0000000000000}*f_curr[13];
+	V x33 = T{72.0000000000000}*f_curr[11];
+	V x44 = T{72.0000000000000}*f_curr[1];
+	V x50 = T{72.0000000000000}*f_curr[17];
+	V x56 = T{72.0000000000000}*f_curr[15];
+	V x64 = T{72.0000000000000}*f_curr[4];
+	V x68 = T{72.0000000000000}*f_curr[14];
+	V x72 = T{72.0000000000000}*f_curr[18];
+	rho = x3;
+	V x78 = T{2.00000000000000}*rho;
+	u[0] = x6*(x0 + x4 + x5 - f_curr[1] - f_curr[5] - f_curr[8]);
+	V x12 = -u[0];
+	V x14 = T{6.00000000000000}*u[0];
+	V x17 = u[0]*u[0];
+	V x18 = T{3.00000000000000}*x17;
+	V x27 = -x14;
+	V x38 = T{2.00000000000000} - x18;
+	V x79 = T{6.00000000000000}*x17;
+	u[1] = x6*(x1 + x4 + x7 - f_curr[12] - f_curr[13] - f_curr[4]);
+	V x10 = T{6.00000000000000}*u[1];
+	V x11 = -x10;
+	V x13 = x12 + u[1];
+	V x19 = u[1]*u[1];
+	V x20 = T{3.00000000000000}*x19;
+	V x21 = x18 + x20 + T{-2.00000000000000};
+	V x28 = -u[1];
+	V x29 = x28 + u[0];
+	V x34 = u[0] + u[1];
+	V x35 = T{9.00000000000000}*(x34*x34);
+	V x37 = -x20;
+	V x39 = x37 + x38;
+	V x85 = T{6.00000000000000}*x19;
+	u[2] = x6*(x2 + x5 + x7 - f_curr[14] - f_curr[16] - f_curr[17]);
+	V x15 = u[2]*u[2];
+	V x16 = T{3.00000000000000}*x15;
+	V x22 = x16 + x21;
+	V x23 = x14 + x22;
+	V x24 = rho*(x11 + x23 - T{9.00000000000000}*x13*x13);
+	V x25 = x24 + x9;
+	V x30 = x10 + x22;
+	V x31 = rho*(x27 + x30 - T{9.00000000000000}*x29*x29);
+	V x32 = x26 + x31;
+	V x36 = rho*(x14 + x30 - x35);
+	V x40 = -x16;
+	V x41 = x10 + x40;
+	V x42 = rho*(x14 + x35 + x39 + x41) - T{72.0000000000000}*f_curr[7];
+	V x43 = -x33 - x36 + x42;
+	V x45 = T{6.00000000000000}*u[2];
+	V x46 = -x45;
+	V x47 = x12 + u[2];
+	V x48 = rho*(x23 + x46 - T{9.00000000000000}*x47*x47);
+	V x49 = x44 + x48;
+	V x51 = -u[2];
+	V x52 = x51 + u[0];
+	V x53 = x22 + x45;
+	V x54 = rho*(x27 + x53 - T{9.00000000000000}*x52*x52);
+	V x55 = x50 + x54;
+	V x57 = u[0] + u[2];
+	V x58 = T{9.00000000000000}*(x57*x57);
+	V x59 = rho*(x14 + x53 - x58);
+	V x60 = x39 + x45;
+	V x61 = x14 + x40;
+	V x62 = rho*(x58 + x60 + x61) - T{72.0000000000000}*f_curr[3];
+	V x63 = -x56 - x59 + x62;
+	V x65 = x28 + u[2];
+	V x66 = rho*(x30 + x46 - T{9.00000000000000}*x65*x65);
+	V x67 = x64 + x66;
+	V x69 = x51 + u[1];
+	V x70 = rho*(x11 + x53 - T{9.00000000000000}*x69*x69);
+	V x71 = x68 + x70;
+	V x73 = u[1] + u[2];
+	V x74 = T{9.00000000000000}*(x73*x73);
+	V x75 = rho*(x30 + x45 - x74);
+	V x76 = rho*(x41 + x60 + x74) - T{72.0000000000000}*f_curr[0];
+	V x77 = -x72 - x75 + x76;
+	V x80 = x16 + T{-2.00000000000000};
+	V x81 = x14 + x20 - x79 + x80;
+	V x82 = x37 + x61 + x79 + T{2.00000000000000};
+	V x83 = -x44 - x48 - x50 - x54 + x63;
+	V x84 = -x24 - x26 - x31 + x43 - x9;
+	V x86 = x10 + x18 + x80 - x85;
+	V x87 = x38 + x41 + x85;
+	V x88 = -x64 - x66 - x68 - x70 + x77;
+	V x89 = T{6.00000000000000}*x15;
+	V x90 = x21 + x45 - x89;
+	V x91 = x60 + x89;
+	V x92 = T{1}/(T{20.0000000000000}*tau + T{2.00000000000000} + T{5.04537849152229}*simd::sqrt<T>((smagorinsky*smagorinsky)*simd::sqrt<T>((x25 + x32 + x43)*(x25 + x32 + x43) + (x49 + x55 + x63)*(x49 + x55 + x63) + (x67 + x71 + x77)*(x67 + x71 + x77) + T{0.500000000000000}*((-x78*x81 + x78*x82 + x83 + x84 - T{72}*f_curr[10] - T{72}*f_curr[8])*(-x78*x81 + x78*x82 + x83 + x84 - T{72}*f_curr[10] - T{72}*f_curr[8])) + T{0.500000000000000}*((-x78*x86 + x78*x87 + x84 + x88 - T{72}*f_curr[12] - T{72}*f_curr[6])*(-x78*x86 + x78*x87 + x84 + x88 - T{72}*f_curr[12] - T{72}*f_curr[6])) + T{0.500000000000000}*((-x78*x90 + x78*x91 + x83 + x88 - T{72}*f_curr[16] - T{72}*f_curr[2])*(-x78*x90 + x78*x91 + x83 + x88 - T{72}*f_curr[16] - T{72}*f_curr[2]))) + T{0.157134840263677}*(x8*x8)));
+	f_next[0] = T{0.333333333333333}*x76*x92 + f_curr[0];
+	f_next[1] = -T{0.333333333333333}*x49*x92 + f_curr[1];
+	f_next[2] = T{0.666666666666667}*x92*(rho*x91 - T{36.0000000000000}*f_curr[2]) + f_curr[2];
+	f_next[3] = T{0.333333333333333}*x62*x92 + f_curr[3];
+	f_next[4] = -T{0.333333333333333}*x67*x92 + f_curr[4];
+	f_next[5] = -T{0.333333333333333}*x25*x92 + f_curr[5];
+	f_next[6] = T{0.666666666666667}*x92*(rho*x87 - T{36.0000000000000}*f_curr[6]) + f_curr[6];
+	f_next[7] = T{0.333333333333333}*x42*x92 + f_curr[7];
+	f_next[8] = -T{0.666666666666667}*x92*(rho*x81 + T{36.0000000000000}*f_curr[8]) + f_curr[8];
+	f_next[9] = -T{4.00000000000000}*x92*(rho*x22 + T{6.00000000000000}*f_curr[9]) + f_curr[9];
+	f_next[10] = T{0.666666666666667}*x92*(rho*x82 - T{36.0000000000000}*f_curr[10]) + f_curr[10];
+	f_next[11] = -T{0.333333333333333}*x92*(x33 + x36) + f_curr[11];
+	f_next[12] = -T{0.666666666666667}*x92*(rho*x86 + T{36.0000000000000}*f_curr[12]) + f_curr[12];
+	f_next[13] = -T{0.333333333333333}*x32*x92 + f_curr[13];
+	f_next[14] = -T{0.333333333333333}*x71*x92 + f_curr[14];
+	f_next[15] = -T{0.333333333333333}*x92*(x56 + x59) + f_curr[15];
+	f_next[16] = -T{0.666666666666667}*x92*(rho*x90 + T{36.0000000000000}*f_curr[16]) + f_curr[16];
+	f_next[17] = -T{0.333333333333333}*x55*x92 + f_curr[17];
+	f_next[18] = -T{0.333333333333333}*x92*(x72 + x75) + f_curr[18];
+}
+
+};
diff --git a/src/LLBM/wall.h b/src/LLBM/wall.h
new file mode 100644
index 0000000..26534e4
--- /dev/null
+++ b/src/LLBM/wall.h
@@ -0,0 +1,4 @@
+#pragma once
+
+template<int N_0, int N_1, int N_2=0>
+struct WallNormal { };