summaryrefslogtreecommitdiff
path: root/tangle
diff options
context:
space:
mode:
Diffstat (limited to 'tangle')
-rw-r--r--tangle/LLBM/base.h6
-rw-r--r--tangle/LLBM/boundary.h13
-rw-r--r--tangle/LLBM/bulk.h3
-rw-r--r--tangle/LLBM/call_tag.h12
-rw-r--r--tangle/LLBM/descriptor.h298
-rw-r--r--tangle/LLBM/kernel/bounce_back.h44
-rw-r--r--tangle/LLBM/kernel/bounce_back_moving_wall.h44
-rw-r--r--tangle/LLBM/kernel/bouzidi.h39
-rw-r--r--tangle/LLBM/kernel/collect_curl.h44
-rw-r--r--tangle/LLBM/kernel/collect_moments.h45
-rw-r--r--tangle/LLBM/kernel/collect_q_criterion.h104
-rw-r--r--tangle/LLBM/kernel/collect_shear_layer_normal.h139
-rw-r--r--tangle/LLBM/kernel/collect_streamlines.h48
-rw-r--r--tangle/LLBM/kernel/collect_velocity_norm.h61
-rw-r--r--tangle/LLBM/kernel/collide.h131
-rw-r--r--tangle/LLBM/kernel/equilibrium_density_wall.h226
-rw-r--r--tangle/LLBM/kernel/equilibrium_velocity_wall.h222
-rw-r--r--tangle/LLBM/kernel/executor.h171
-rw-r--r--tangle/LLBM/kernel/free_slip.h82
-rw-r--r--tangle/LLBM/kernel/initialize.h44
-rw-r--r--tangle/LLBM/kernel/propagate.h20
-rw-r--r--tangle/LLBM/kernel/smagorinsky_collide.h183
-rw-r--r--tangle/LLBM/lattice.h131
-rw-r--r--tangle/LLBM/materials.h108
-rw-r--r--tangle/LLBM/memory.h134
-rw-r--r--tangle/LLBM/operator.h27
-rw-r--r--tangle/LLBM/propagate.h111
-rw-r--r--tangle/LLBM/sdf.h75
-rw-r--r--tangle/LLBM/sdf_boundary.h114
-rw-r--r--tangle/LLBM/volumetric.h134
-rw-r--r--tangle/LLBM/wall.h4
-rw-r--r--tangle/asset/noise/blue_0.pngbin0 -> 3763 bytes
-rw-r--r--tangle/asset/noise/blue_1.pngbin0 -> 3775 bytes
-rw-r--r--tangle/asset/noise/blue_2.pngbin0 -> 3752 bytes
-rw-r--r--tangle/asset/noise/blue_3.pngbin0 -> 3789 bytes
-rw-r--r--tangle/asset/noise/blue_4.pngbin0 -> 3754 bytes
-rw-r--r--tangle/asset/palette/4wave_ROTB.pngbin0 -> 1572 bytes
-rw-r--r--tangle/asset/palette/4wave_equal.pngbin0 -> 1485 bytes
-rw-r--r--tangle/asset/palette/5wave_cool.pngbin0 -> 1601 bytes
-rw-r--r--tangle/asset/palette/autumn.pngbin0 -> 827 bytes
-rw-r--r--tangle/asset/palette/blue.pngbin0 -> 1223 bytes
-rw-r--r--tangle/asset/palette/blue_orange.pngbin0 -> 1457 bytes
-rw-r--r--tangle/asset/palette/green_brown.pngbin0 -> 1269 bytes
-rw-r--r--tangle/asset/palette/orange.pngbin0 -> 1165 bytes
-rw-r--r--tangle/asset/shader/blur.frag20
-rw-r--r--tangle/benchmark-ldc.cu95
-rw-r--r--tangle/channel-with-sphere.cu85
-rw-r--r--tangle/ldc-2d.cu82
-rw-r--r--tangle/ldc-3d.cu58
-rw-r--r--tangle/magnus.cu116
-rw-r--r--tangle/nozzle.cu73
-rw-r--r--tangle/sampler/curl_norm.h77
-rw-r--r--tangle/sampler/q_criterion.h90
-rw-r--r--tangle/sampler/sampler.h32
-rw-r--r--tangle/sampler/shear_layer.h67
-rw-r--r--tangle/sampler/velocity_norm.h79
-rw-r--r--tangle/taylor-couette.cu75
-rw-r--r--tangle/tmp/noise_overview.pngbin0 -> 14996 bytes
-rw-r--r--tangle/tmp/test_noise.pngbin0 -> 26571 bytes
-rw-r--r--tangle/util/camera.h86
-rw-r--r--tangle/util/colormap.h38
-rw-r--r--tangle/util/noise.h39
-rw-r--r--tangle/util/render_window.h95
-rw-r--r--tangle/util/texture.h25
-rw-r--r--tangle/util/timer.h19
-rw-r--r--tangle/util/volumetric_example.h121
66 files changed, 4189 insertions, 0 deletions
diff --git a/tangle/LLBM/base.h b/tangle/LLBM/base.h
new file mode 100644
index 0000000..e356e02
--- /dev/null
+++ b/tangle/LLBM/base.h
@@ -0,0 +1,6 @@
+#pragma once
+
+#include "descriptor.h"
+#include "memory.h"
+#include "lattice.h"
+#include "materials.h"
diff --git a/tangle/LLBM/boundary.h b/tangle/LLBM/boundary.h
new file mode 100644
index 0000000..08793a0
--- /dev/null
+++ b/tangle/LLBM/boundary.h
@@ -0,0 +1,13 @@
+#include "kernel/bounce_back.h"
+
+#include "kernel/bounce_back_moving_wall.h"
+
+#include "kernel/free_slip.h"
+
+#include "kernel/bouzidi.h"
+
+#include "kernel/equilibrium_velocity_wall.h"
+
+#include "kernel/equilibrium_density_wall.h"
+
+#include "sdf_boundary.h"
diff --git a/tangle/LLBM/bulk.h b/tangle/LLBM/bulk.h
new file mode 100644
index 0000000..62d704a
--- /dev/null
+++ b/tangle/LLBM/bulk.h
@@ -0,0 +1,3 @@
+#include "kernel/collide.h"
+
+#include "kernel/smagorinsky_collide.h"
diff --git a/tangle/LLBM/call_tag.h b/tangle/LLBM/call_tag.h
new file mode 100644
index 0000000..fe559ef
--- /dev/null
+++ b/tangle/LLBM/call_tag.h
@@ -0,0 +1,12 @@
+#pragma once
+
+namespace tag {
+
+struct call_by_cell_id { };
+struct call_by_list_index { };
+struct call_by_spatial_cell_mask { };
+
+struct post_process_by_list_index { };
+struct post_process_by_spatial_cell_mask { };
+
+}
diff --git a/tangle/LLBM/descriptor.h b/tangle/LLBM/descriptor.h
new file mode 100644
index 0000000..d4ca44d
--- /dev/null
+++ b/tangle/LLBM/descriptor.h
@@ -0,0 +1,298 @@
+#pragma once
+
+#include <algorithm>
+#include <cstdint>
+#include <type_traits>
+#include <cuda-samples/Common/helper_math.h>
+
+#ifdef __CUDA_ARCH__
+ #define DATA device_data
+#else
+ #define DATA host_data
+#endif
+
+using pop_index_t = std::uint8_t;
+
+namespace descriptor {
+
+struct D2Q9 {
+ static constexpr unsigned d = 2;
+ static constexpr unsigned q = 9;
+};
+
+struct D3Q19 {
+ static constexpr unsigned d = 3;
+ static constexpr unsigned q = 19;
+};
+
+namespace device_data {
+ template <typename DESCRIPTOR>
+ __constant__ pop_index_t opposite[DESCRIPTOR::q] { };
+
+ template <typename DESCRIPTOR>
+ __constant__ int c[DESCRIPTOR::q][DESCRIPTOR::d] { };
+
+ template <typename DESCRIPTOR>
+ __constant__ float c_length[DESCRIPTOR::q] { };
+
+ template <typename DESCRIPTOR>
+ __constant__ float weight[DESCRIPTOR::q] { };
+
+ template <>
+ __constant__ pop_index_t opposite<D2Q9>[9] = {
+ 8, 7, 6, 5, 4, 3, 2, 1, 0
+ };
+
+ template <>
+ __constant__ int c<D2Q9>[9][2] = {
+ {-1,-1}, {-1,0}, {-1,1}, {0,-1}, {0,0}, {0,1}, {1,-1}, {1,0}, {1,1}
+ };
+
+ template <>
+ __constant__ float c_length<D2Q9>[9] = {
+ 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.00000000000000, 0, 1.00000000000000, 1.41421356237310, 1.00000000000000, 1.41421356237310
+ };
+
+ template <>
+ __constant__ float weight<D2Q9>[9] = {
+ 0.0277777777777778, 0.111111111111111, 0.0277777777777778, 0.111111111111111, 0.444444444444444, 0.111111111111111, 0.0277777777777778, 0.111111111111111, 0.0277777777777778
+ };
+
+ template <>
+ __constant__ pop_index_t opposite<D3Q19>[19] = {
+ 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0
+ };
+
+ template <>
+ __constant__ int c<D3Q19>[19][3] = {
+ {0,1,1}, {-1,0,1}, {0,0,1}, {1,0,1}, {0,-1,1}, {-1,1,0}, {0,1,0}, {1,1,0}, {-1,0,0}, {0,0,0}, {1,0,0}, {-1,-1,0}, {0,-1,0}, {1,-1,0}, {0,1,-1}, {-1,0,-1}, {0,0,-1}, {1,0,-1}, {0,-1,-1}
+ };
+
+ template <>
+ __constant__ float c_length<D3Q19>[19] = {
+ 1.41421356237310, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.41421356237310, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.00000000000000, 0, 1.00000000000000, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.41421356237310, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.41421356237310
+ };
+
+ template <>
+ __constant__ float weight<D3Q19>[19] = {
+ 0.0277777777777778, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0277777777777778, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0555555555555556, 0.333333333333333, 0.0555555555555556, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0277777777777778, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0277777777777778
+ };
+}
+
+namespace host_data {
+ template <typename DESCRIPTOR>
+ constexpr pop_index_t opposite[DESCRIPTOR::q] { };
+
+ template <typename DESCRIPTOR>
+ constexpr int c[DESCRIPTOR::q][DESCRIPTOR::d] { };
+
+ template <typename DESCRIPTOR>
+ constexpr float c_length[DESCRIPTOR::q] { };
+
+ template <typename DESCRIPTOR>
+ constexpr float weight[DESCRIPTOR::q] { };
+
+ template <>
+ constexpr pop_index_t opposite<D2Q9>[9] = {
+ 8, 7, 6, 5, 4, 3, 2, 1, 0
+ };
+
+ template <>
+ constexpr int c<D2Q9>[9][2] = {
+ {-1,-1}, {-1,0}, {-1,1}, {0,-1}, {0,0}, {0,1}, {1,-1}, {1,0}, {1,1}
+ };
+
+ template <>
+ constexpr float c_length<D2Q9>[9] = {
+ 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.00000000000000, 0, 1.00000000000000, 1.41421356237310, 1.00000000000000, 1.41421356237310
+ };
+
+ template <>
+ constexpr float weight<D2Q9>[9] = {
+ 0.0277777777777778, 0.111111111111111, 0.0277777777777778, 0.111111111111111, 0.444444444444444, 0.111111111111111, 0.0277777777777778, 0.111111111111111, 0.0277777777777778
+ };
+
+ template <>
+ constexpr pop_index_t opposite<D3Q19>[19] = {
+ 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0
+ };
+
+ template <>
+ constexpr int c<D3Q19>[19][3] = {
+ {0,1,1}, {-1,0,1}, {0,0,1}, {1,0,1}, {0,-1,1}, {-1,1,0}, {0,1,0}, {1,1,0}, {-1,0,0}, {0,0,0}, {1,0,0}, {-1,-1,0}, {0,-1,0}, {1,-1,0}, {0,1,-1}, {-1,0,-1}, {0,0,-1}, {1,0,-1}, {0,-1,-1}
+ };
+
+ template <>
+ constexpr float c_length<D3Q19>[19] = {
+ 1.41421356237310, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.41421356237310, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.00000000000000, 0, 1.00000000000000, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.41421356237310, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.41421356237310
+ };
+
+ template <>
+ constexpr float weight<D3Q19>[19] = {
+ 0.0277777777777778, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0277777777777778, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0555555555555556, 0.333333333333333, 0.0555555555555556, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0277777777777778, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0277777777777778
+ };
+}
+
+template <typename DESCRIPTOR>
+__host__ __device__
+pop_index_t opposite(pop_index_t iPop) {
+ return DESCRIPTOR::q - 1 - iPop;
+}
+
+template <typename DESCRIPTOR>
+__host__ __device__
+int velocity(pop_index_t iPop, unsigned iDim) {
+ return DATA::template c<DESCRIPTOR>[iPop][iDim];
+}
+
+template <typename DESCRIPTOR>
+__host__ __device__
+std::enable_if_t<DESCRIPTOR::d == 2, float2> velocity(pop_index_t iPop) {
+ return make_float2(DATA::template c<DESCRIPTOR>[iPop][0],
+ DATA::template c<DESCRIPTOR>[iPop][1]);
+}
+
+template <typename DESCRIPTOR>
+__host__ __device__
+std::enable_if_t<DESCRIPTOR::d == 3, float3> velocity(pop_index_t iPop) {
+ return make_float3(DATA::template c<DESCRIPTOR>[iPop][0],
+ DATA::template c<DESCRIPTOR>[iPop][1],
+ DATA::template c<DESCRIPTOR>[iPop][2]);
+}
+
+template <typename DESCRIPTOR>
+__host__ __device__
+float velocity_length(pop_index_t iPop) {
+ return DATA::template c_length<DESCRIPTOR>[iPop];
+}
+
+template <typename DESCRIPTOR>
+__host__ __device__
+float weight(pop_index_t iPop) {
+ return DATA::template weight<DESCRIPTOR>[iPop];
+}
+
+template <unsigned D>
+struct CuboidD;
+
+template <>
+struct CuboidD<2> {
+ const std::size_t nX;
+ const std::size_t nY;
+ const std::size_t nZ;
+ const std::size_t volume;
+
+ CuboidD(std::size_t x, std::size_t y):
+ nX(x), nY(y), nZ(1),
+ volume(x*y) { };
+};
+
+template <>
+struct CuboidD<3> {
+ const std::size_t nX;
+ const std::size_t nY;
+ const std::size_t nZ;
+ const std::size_t volume;
+ const std::size_t plane;
+
+ CuboidD(std::size_t x, std::size_t y, std::size_t z):
+ nX(x), nY(y), nZ(z),
+ volume(x*y*z),
+ plane(x*y) { };
+};
+
+template <typename DESCRIPTOR>
+using Cuboid = CuboidD<DESCRIPTOR::d>;
+
+__host__ __device__
+std::size_t gid(const CuboidD<2>& c, int iX, int iY, int iZ=0) {
+ return iY*c.nX + iX;
+}
+
+__host__ __device__
+std::size_t gid(const CuboidD<3>& c, int iX, int iY, int iZ) {
+ return iZ*c.plane + iY*c.nX + iX;
+}
+
+__host__ __device__
+int offset(const CuboidD<2>& c, int iX, int iY) {
+ return iY*c.nX + iX;
+}
+
+template <typename DESCRIPTOR>
+__host__ __device__
+int offset(const CuboidD<2>& c, pop_index_t iPop) {
+ static_assert(DESCRIPTOR::d == 2, "Dimensions must match");
+ return offset(c,
+ descriptor::velocity<DESCRIPTOR>(iPop, 0),
+ descriptor::velocity<DESCRIPTOR>(iPop, 1)
+ );
+}
+
+__host__ __device__
+int offset(const CuboidD<3>& c, int iX, int iY, int iZ) {
+ return iZ*c.plane + iY*c.nX + iX;
+}
+
+template <typename DESCRIPTOR>
+__host__ __device__
+int offset(const CuboidD<3>& c, pop_index_t iPop) {
+ static_assert(DESCRIPTOR::d == 3, "Dimensions must match");
+ return offset(c,
+ descriptor::velocity<DESCRIPTOR>(iPop, 0),
+ descriptor::velocity<DESCRIPTOR>(iPop, 1),
+ descriptor::velocity<DESCRIPTOR>(iPop, 2)
+ );
+}
+
+template <typename DESCRIPTOR>
+__host__ __device__
+std::size_t neighbor(const CuboidD<2>& c, std::size_t iCell, pop_index_t iPop) {
+ return iCell + offset<DESCRIPTOR>(c, iPop);
+}
+
+template <typename DESCRIPTOR>
+__host__ __device__
+std::size_t neighbor(const CuboidD<3>& c, std::size_t iCell, pop_index_t iPop) {
+ return iCell + offset<DESCRIPTOR>(c, iPop);
+}
+
+__host__ __device__
+uint2 gidInverse(const CuboidD<2>& c, std::size_t gid) {
+ int iY = gid / c.nX;
+ int iX = gid % c.nX;
+ return make_uint2(iX, iY);
+}
+
+__host__ __device__
+uint3 gidInverse(const CuboidD<3>& c, std::size_t gid) {
+ int iZ = gid / c.plane;
+ int iY = (gid % c.plane) / c.nX;
+ int iX = (gid % c.plane) % c.nX;
+ return make_uint3(iX,iY,iZ);
+}
+
+__host__ __device__
+float2 gidInverseSmooth(const CuboidD<2>& c, std::size_t gid) {
+ int iY = gid / c.nX;
+ int iX = gid % c.nX;
+ return make_float2(iX, iY);
+}
+
+__host__ __device__
+float3 gidInverseSmooth(const CuboidD<3>& c, std::size_t gid) {
+ int iZ = gid / c.plane;
+ int iY = (gid % c.plane) / c.nX;
+ int iX = (gid % c.plane) % c.nX;
+ return make_float3(iX,iY,iZ);
+}
+
+bool isInside(const CuboidD<2>& c, std::size_t gid) {
+ return gid < c.volume;
+}
+
+bool isInside(const CuboidD<3>& c, std::size_t gid) {
+ return gid < c.volume;
+}
+
+}
diff --git a/tangle/LLBM/kernel/bounce_back.h b/tangle/LLBM/kernel/bounce_back.h
new file mode 100644
index 0000000..d7ec4a7
--- /dev/null
+++ b/tangle/LLBM/kernel/bounce_back.h
@@ -0,0 +1,44 @@
+#pragma once
+#include <LLBM/call_tag.h>
+
+struct BounceBackO {
+
+using call_tag = tag::call_by_cell_id;
+
+template <typename T, typename S>
+__device__ static void apply(descriptor::D2Q9, S f_curr[9], S f_next[9], std::size_t gid) {
+ f_next[0] = f_curr[8];
+ f_next[1] = f_curr[7];
+ f_next[2] = f_curr[6];
+ f_next[3] = f_curr[5];
+ f_next[4] = f_curr[4];
+ f_next[5] = f_curr[3];
+ f_next[6] = f_curr[2];
+ f_next[7] = f_curr[1];
+ f_next[8] = f_curr[0];
+}
+
+template <typename T, typename S>
+__device__ static void apply(descriptor::D3Q19, S f_curr[19], S f_next[19], std::size_t gid) {
+ 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/tangle/LLBM/kernel/bounce_back_moving_wall.h b/tangle/LLBM/kernel/bounce_back_moving_wall.h
new file mode 100644
index 0000000..c4b40f9
--- /dev/null
+++ b/tangle/LLBM/kernel/bounce_back_moving_wall.h
@@ -0,0 +1,44 @@
+#pragma once
+#include <LLBM/call_tag.h>
+
+struct BounceBackMovingWallO {
+
+using call_tag = tag::call_by_cell_id;
+
+template <typename T, typename S>
+__device__ static void apply(descriptor::D2Q9, S f_curr[9], S f_next[9], std::size_t gid, T u_0, T u_1) {
+ f_next[0] = -T{0.166666666666667}*u_0 - T{0.166666666666667}*u_1 + f_curr[8];
+ f_next[1] = -T{0.666666666666667}*u_0 + f_curr[7];
+ f_next[2] = -T{0.166666666666667}*u_0 + T{0.166666666666667}*u_1 + f_curr[6];
+ f_next[3] = -T{0.666666666666667}*u_1 + f_curr[5];
+ f_next[4] = f_curr[4];
+ f_next[5] = T{0.666666666666667}*u_1 + f_curr[3];
+ f_next[6] = T{0.166666666666667}*u_0 - T{0.166666666666667}*u_1 + f_curr[2];
+ f_next[7] = T{0.666666666666667}*u_0 + f_curr[1];
+ f_next[8] = T{0.166666666666667}*u_0 + T{0.166666666666667}*u_1 + f_curr[0];
+}
+
+template <typename T, typename S>
+__device__ static void apply(descriptor::D3Q19, S f_curr[19], S f_next[19], std::size_t gid, T u_0, T u_1, T u_2) {
+ f_next[0] = T{0.166666666666667}*u_1 + T{0.166666666666667}*u_2 + f_curr[18];
+ f_next[1] = -T{0.166666666666667}*u_0 + T{0.166666666666667}*u_2 + f_curr[17];
+ f_next[2] = T{0.333333333333333}*u_2 + f_curr[16];
+ f_next[3] = T{0.166666666666667}*u_0 + T{0.166666666666667}*u_2 + f_curr[15];
+ f_next[4] = -T{0.166666666666667}*u_1 + T{0.166666666666667}*u_2 + f_curr[14];
+ f_next[5] = -T{0.166666666666667}*u_0 + T{0.166666666666667}*u_1 + f_curr[13];
+ f_next[6] = T{0.333333333333333}*u_1 + f_curr[12];
+ f_next[7] = T{0.166666666666667}*u_0 + T{0.166666666666667}*u_1 + f_curr[11];
+ f_next[8] = -T{0.333333333333333}*u_0 + f_curr[10];
+ f_next[9] = f_curr[9];
+ f_next[10] = T{0.333333333333333}*u_0 + f_curr[8];
+ f_next[11] = -T{0.166666666666667}*u_0 - T{0.166666666666667}*u_1 + f_curr[7];
+ f_next[12] = -T{0.333333333333333}*u_1 + f_curr[6];
+ f_next[13] = T{0.166666666666667}*u_0 - T{0.166666666666667}*u_1 + f_curr[5];
+ f_next[14] = T{0.166666666666667}*u_1 - T{0.166666666666667}*u_2 + f_curr[4];
+ f_next[15] = -T{0.166666666666667}*u_0 - T{0.166666666666667}*u_2 + f_curr[3];
+ f_next[16] = -T{0.333333333333333}*u_2 + f_curr[2];
+ f_next[17] = T{0.166666666666667}*u_0 - T{0.166666666666667}*u_2 + f_curr[1];
+ f_next[18] = -T{0.166666666666667}*u_1 - T{0.166666666666667}*u_2 + f_curr[0];
+}
+
+};
diff --git a/tangle/LLBM/kernel/bouzidi.h b/tangle/LLBM/kernel/bouzidi.h
new file mode 100644
index 0000000..b1abd2c
--- /dev/null
+++ b/tangle/LLBM/kernel/bouzidi.h
@@ -0,0 +1,39 @@
+#pragma once
+#include <LLBM/call_tag.h>
+#include <LLBM/lattice.h>
+
+template <typename S>
+struct BouzidiConfig {
+ std::size_t* boundary; // boundary cell to be interpolated
+ std::size_t* solid; // adjacent solid cell
+ std::size_t* fluid; // adjacent fluid cell
+ S* distance; // precomputed distance factor q
+ S* correction; // correction for moving walls
+ pop_index_t* missing; // population to be reconstructed
+};
+
+struct BouzidiO {
+
+using call_tag = tag::call_by_list_index;
+
+template <typename T, typename S, typename DESCRIPTOR>
+__device__ static void apply(
+ LatticeView<DESCRIPTOR,S> lattice
+ , std::size_t index
+ , std::size_t count
+ , BouzidiConfig<S> config
+) {
+ pop_index_t& iPop = config.missing[index];
+ pop_index_t jPop = descriptor::opposite<DESCRIPTOR>(iPop);
+ pop_index_t kPop = config.boundary[index] == config.fluid[index] ? iPop : jPop;
+
+ S f_bound_j = *lattice.pop(jPop, config.boundary[index]);
+ S f_fluid_j = *lattice.pop(kPop, config.fluid[index]);
+ S* f_next_i = lattice.pop(iPop, config.solid[index]);
+
+ *f_next_i = config.distance[index] * f_bound_j
+ + (1. - config.distance[index]) * f_fluid_j
+ + config.correction[index];
+}
+
+};
diff --git a/tangle/LLBM/kernel/collect_curl.h b/tangle/LLBM/kernel/collect_curl.h
new file mode 100644
index 0000000..054b359
--- /dev/null
+++ b/tangle/LLBM/kernel/collect_curl.h
@@ -0,0 +1,44 @@
+#pragma once
+#include <LLBM/call_tag.h>
+
+struct CollectCurlF {
+
+using call_tag = tag::call_by_spatial_cell_mask;
+
+template <typename T, typename S>
+__device__ static void apply(
+ descriptor::D3Q19
+ , S f_curr[19]
+ , descriptor::CuboidD<3> cuboid
+ , std::size_t gid
+ , std::size_t iX
+ , std::size_t iY
+ , std::size_t iZ
+ , S* moments_u
+ , cudaSurfaceObject_t surface
+ , S* curl_norm = nullptr
+) {
+ auto u_x = [moments_u,cuboid,gid] __device__ (int x, int y, int z) -> T {
+ return moments_u[3*(gid + descriptor::offset(cuboid,x,y,z)) + 0];
+ };
+ auto u_y = [moments_u,cuboid,gid] __device__ (int x, int y, int z) -> T {
+ return moments_u[3*(gid + descriptor::offset(cuboid,x,y,z)) + 1];
+ };
+ auto u_z = [moments_u,cuboid,gid] __device__ (int x, int y, int z) -> T {
+ return moments_u[3*(gid + descriptor::offset(cuboid,x,y,z)) + 2];
+ };
+
+ T curl_0 = T{0.500000000000000}*u_y(0, 0, -1) - T{0.500000000000000}*u_y(0, 0, 1) - T{0.500000000000000}*u_z(0, -1, 0) + T{0.500000000000000}*u_z(0, 1, 0);
+ T curl_1 = -T{0.500000000000000}*u_x(0, 0, -1) + T{0.500000000000000}*u_x(0, 0, 1) + T{0.500000000000000}*u_z(-1, 0, 0) - T{0.500000000000000}*u_z(1, 0, 0);
+ T curl_2 = T{0.500000000000000}*u_x(0, -1, 0) - T{0.500000000000000}*u_x(0, 1, 0) - T{0.500000000000000}*u_y(-1, 0, 0) + T{0.500000000000000}*u_y(1, 0, 0);
+ float3 curl = make_float3(curl_0, curl_1, curl_2);
+ float norm = length(curl);
+
+ surf3Dwrite(norm, surface, iX*sizeof(float), iY, iZ);
+
+ if (curl_norm != nullptr) {
+ curl_norm[gid] = norm;
+ }
+}
+
+};
diff --git a/tangle/LLBM/kernel/collect_moments.h b/tangle/LLBM/kernel/collect_moments.h
new file mode 100644
index 0000000..b2ef09b
--- /dev/null
+++ b/tangle/LLBM/kernel/collect_moments.h
@@ -0,0 +1,45 @@
+#pragma once
+#include <LLBM/call_tag.h>
+
+struct CollectMomentsF {
+
+using call_tag = tag::call_by_cell_id;
+
+template <typename T, typename S>
+__device__ static void apply(descriptor::D2Q9, S f_curr[9], std::size_t gid, T* cell_rho, T* cell_u) {
+ T m0 = f_curr[1] + f_curr[2];
+ T m1 = f_curr[3] + f_curr[6];
+ T m2 = m0 + m1 + f_curr[0] + f_curr[4] + f_curr[5] + f_curr[7] + f_curr[8];
+ T m3 = f_curr[0] - f_curr[8];
+ T m4 = T{1} / (m2);
+ T rho = m2;
+ T u_0 = -m4*(m0 + m3 - f_curr[6] - f_curr[7]);
+ T u_1 = -m4*(m1 + m3 - f_curr[2] - f_curr[5]);
+
+ cell_rho[gid] = rho;
+ cell_u[2*gid+0] = u_0;
+ cell_u[2*gid+1] = u_1;
+}
+
+template <typename T, typename S>
+__device__ static void apply(descriptor::D3Q19, S f_curr[19], std::size_t gid, T* cell_rho, T* cell_u) {
+ T m0 = f_curr[10] + f_curr[13] + f_curr[17];
+ T m1 = f_curr[14] + f_curr[5] + f_curr[6];
+ T m2 = f_curr[1] + f_curr[2] + f_curr[4];
+ T 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];
+ T m4 = -f_curr[11] + f_curr[7];
+ T m5 = -f_curr[15] + f_curr[3];
+ T m6 = T{1} / (m3);
+ T m7 = f_curr[0] - f_curr[18];
+ T rho = m3;
+ T u_0 = m6*(m0 + m4 + m5 - f_curr[1] - f_curr[5] - f_curr[8]);
+ T u_1 = m6*(m1 + m4 + m7 - f_curr[12] - f_curr[13] - f_curr[4]);
+ T u_2 = m6*(m2 + m5 + m7 - f_curr[14] - f_curr[16] - f_curr[17]);
+
+ cell_rho[gid] = rho;
+ cell_u[3*gid+0] = u_0;
+ cell_u[3*gid+1] = u_1;
+ cell_u[3*gid+2] = u_2;
+}
+
+};
diff --git a/tangle/LLBM/kernel/collect_q_criterion.h b/tangle/LLBM/kernel/collect_q_criterion.h
new file mode 100644
index 0000000..19b7f68
--- /dev/null
+++ b/tangle/LLBM/kernel/collect_q_criterion.h
@@ -0,0 +1,104 @@
+#pragma once
+#include <LLBM/call_tag.h>
+
+struct CollectQCriterionF {
+
+using call_tag = tag::call_by_spatial_cell_mask;
+
+template <typename T, typename S>
+__device__ static void apply(
+ descriptor::D3Q19
+ , S f_curr[19]
+ , descriptor::CuboidD<3> cuboid
+ , std::size_t gid
+ , std::size_t iX
+ , std::size_t iY
+ , std::size_t iZ
+ , T* cell_rho
+ , T* cell_u
+ , T* cell_curl_norm
+ , cudaSurfaceObject_t surface
+ , T* cell_q = nullptr
+) {
+ const T rho = cell_rho[gid];
+ const T u_0 = cell_u[3*gid + 0];
+ const T u_1 = cell_u[3*gid + 1];
+ const T u_2 = cell_u[3*gid + 2];
+
+ T x0 = T{72.0000000000000}*f_curr[5];
+ T x1 = T{72.0000000000000}*f_curr[13];
+ T x2 = T{6.00000000000000}*u_1;
+ T x3 = -x2;
+ T x4 = -u_0;
+ T x5 = x4 + u_1;
+ T x6 = u_1*u_1;
+ T x7 = T{3.00000000000000}*x6;
+ T x8 = u_0*u_0;