Extract first public LiterateLB version

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;
+  T x9 = T{3.00000000000000}*x8;
+  T x10 = x9 + T{-2.00000000000000};
+  T x11 = x10 + x7;
+  T x12 = u_2*u_2;
+  T x13 = T{3.00000000000000}*x12;
+  T x14 = T{6.00000000000000}*u_0;
+  T x15 = x13 + x14;
+  T x16 = rho*(x11 + x15 + x3 - T{9.00000000000000}*x5*x5);
+  T x17 = -x14;
+  T x18 = -u_1;
+  T x19 = x18 + u_0;
+  T x20 = x13 + x2;
+  T x21 = x11 + x20;
+  T x22 = rho*(x17 + x21 - T{9.00000000000000}*x19*x19);
+  T x23 = u_0 + u_1;
+  T x24 = T{9.00000000000000}*(x23*x23);
+  T x25 = -x7;
+  T x26 = T{2.00000000000000} - x9;
+  T x27 = x25 + x26;
+  T x28 = -x13;
+  T x29 = x2 + x28;
+  T x30 = -rho*(x14 + x21 - x24) + rho*(x14 + x24 + x27 + x29) - T{72.0000000000000}*f_curr[11] - T{72.0000000000000}*f_curr[7];
+  T x31 = T{72.0000000000000}*f_curr[1];
+  T x32 = T{72.0000000000000}*f_curr[17];
+  T x33 = x4 + u_2;
+  T x34 = T{6.00000000000000}*u_2;
+  T x35 = x11 - x34;
+  T x36 = rho*(x15 + x35 - T{9.00000000000000}*x33*x33);
+  T x37 = -u_2;
+  T x38 = x37 + u_0;
+  T x39 = x11 + x34;
+  T x40 = x13 + x39;
+  T x41 = rho*(x17 + x40 - T{9.00000000000000}*x38*x38);
+  T x42 = u_0 + u_2;
+  T x43 = T{9.00000000000000}*(x42*x42);
+  T x44 = x27 + x34;
+  T x45 = x14 + x28;
+  T x46 = -rho*(x15 + x39 - x43) + rho*(x43 + x44 + x45) - T{72.0000000000000}*f_curr[15] - T{72.0000000000000}*f_curr[3];
+  T x47 = T{72.0000000000000}*f_curr[4];
+  T x48 = T{72.0000000000000}*f_curr[14];
+  T x49 = x18 + u_2;
+  T x50 = rho*(x20 + x35 - T{9.00000000000000}*x49*x49);
+  T x51 = x37 + u_1;
+  T x52 = rho*(x3 + x40 - T{9.00000000000000}*x51*x51);
+  T x53 = u_1 + u_2;
+  T x54 = T{9.00000000000000}*(x53*x53);
+  T x55 = -rho*(x20 + x39 - x54) + rho*(x29 + x44 + x54) - T{72.0000000000000}*f_curr[0] - T{72.0000000000000}*f_curr[18];
+  T x56 = T{2.00000000000000}*rho;
+  T x57 = T{6.00000000000000}*x8;
+  T x58 = -x31 - x32 - x36 - x41 + x46;
+  T x59 = -x0 - x1 - x16 - x22 + x30;
+  T x60 = T{6.00000000000000}*x6;
+  T x61 = -x47 - x48 - x50 - x52 + x55;
+  T x62 = T{6.00000000000000}*x12;
+  T strain = T{0.0277777777777778}*sqrt((x0 + x1 + x16 + x22 + x30)*(x0 + x1 + x16 + x22 + x30) + (x31 + x32 + x36 + x41 + x46)*(x31 + x32 + x36 + x41 + x46) + (x47 + x48 + x50 + x52 + x55)*(x47 + x48 + x50 + x52 + x55) + T{0.500000000000000}*((-x56*(x39 - x62) + x56*(x44 + x62) + x58 + x61 - 72*f_curr[16] - 72*f_curr[2])*(-x56*(x39 - x62) + x56*(x44 + x62) + x58 + x61 - 72*f_curr[16] - 72*f_curr[2])) + T{0.500000000000000}*((-x56*(x10 + x20 - x60) + x56*(x26 + x29 + x60) + x59 + x61 - 72*f_curr[12] - 72*f_curr[6])*(-x56*(x10 + x20 - x60) + x56*(x26 + x29 + x60) + x59 + x61 - 72*f_curr[12] - 72*f_curr[6])) + T{0.500000000000000}*((-x56*(x15 - x57 + x7 - 2) + x56*(x25 + x45 + x57 + 2) + x58 + x59 - 72*f_curr[10] - 72*f_curr[8])*(-x56*(x15 - x57 + x7 - 2) + x56*(x25 + x45 + x57 + 2) + x58 + x59 - 72*f_curr[10] - 72*f_curr[8])));
+
+  float vorticity = cell_curl_norm[gid];
+  float q = vorticity*vorticity - strain*strain;
+  q = q > 0 ? q : 0;
+
+  surf3Dwrite(q, surface, iX*sizeof(float), iY, iZ);
+
+  if (cell_q != nullptr) {
+    cell_q[gid] = q;
+  }
+}
+
+};
diff --git a/tangle/LLBM/kernel/collect_shear_layer_normal.h b/tangle/LLBM/kernel/collect_shear_layer_normal.h
new file mode 100644
index 0000000..7bf6eff
--- /dev/null
+++ b/tangle/LLBM/kernel/collect_shear_layer_normal.h
@@ -0,0 +1,139 @@
+#pragma once
+#include <LLBM/call_tag.h>
+
+struct CollectShearLayerNormalsF {
+
+using call_tag = tag::call_by_cell_id;
+
+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* cell_shear_normal
+) {
+  T x0 = f_curr[10] + f_curr[13] + f_curr[17];
+  T x1 = f_curr[14] + f_curr[5] + f_curr[6];
+  T x2 = f_curr[1] + f_curr[2] + f_curr[4];
+  T 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];
+  T x4 = -f_curr[11] + f_curr[7];
+  T x5 = -f_curr[15] + f_curr[3];
+  T x6 = T{1} / (x3);
+  T x7 = f_curr[0] - f_curr[18];
+  T x14 = T{72.0000000000000}*f_curr[5];
+  T x15 = T{72.0000000000000}*f_curr[13];
+  T x39 = T{72.0000000000000}*f_curr[1];
+  T x40 = T{72.0000000000000}*f_curr[17];
+  T x61 = T{72.0000000000000}*f_curr[4];
+  T x62 = T{72.0000000000000}*f_curr[14];
+  T rho = x3;
+  T x56 = T{2.00000000000000}*rho;
+  T u_0 = x6*(x0 + x4 + x5 - f_curr[1] - f_curr[5] - f_curr[8]);
+  T x9 = u_0*u_0;
+  T x16 = -u_0;
+  T x18 = -T{3.00000000000000}*x9;
+  T x21 = T{6.00000000000000}*u_0;
+  T x22 = -x21;
+  T x55 = T{0.0277777777777778}*u_0;
+  T u_1 = x6*(x1 + x4 + x7 - f_curr[12] - f_curr[13] - f_curr[4]);
+  T x8 = T{0.0277777777777778}*u_1;
+  T x10 = u_1*u_1;
+  T x17 = x16 + u_1;
+  T x19 = T{6.00000000000000}*u_1;
+  T x20 = x18 + x19;
+  T x23 = -T{3.00000000000000}*x10;
+  T x28 = -u_1;
+  T x29 = x28 + u_0;
+  T x30 = x18 - x19;
+  T x33 = u_0 + u_1;
+  T x34 = T{9.00000000000000}*(x33*x33);
+  T u_2 = x6*(x2 + x5 + x7 - f_curr[14] - f_curr[16] - f_curr[17]);
+  T x11 = u_2*u_2;
+  T x12 = x10 + x11 + x9;
+  T x13 = pow(x12, T{-0.500000000000000});
+  T x24 = T{2.00000000000000} - T{3.00000000000000}*x11;
+  T x25 = x23 + x24;
+  T x26 = x22 + x25;
+  T x27 = rho*(x20 + x26 + T{9.00000000000000}*(x17*x17));
+  T x31 = x21 + x25;
+  T x32 = rho*(x30 + x31 + T{9.00000000000000}*(x29*x29));
+  T x35 = rho*(x20 + x31 + x34) + rho*(x26 + x30 + x34) - T{72.0000000000000}*f_curr[11] - T{72.0000000000000}*f_curr[7];
+  T x36 = x14 + x15 - x27 - x32 + x35;
+  T x37 = x13*x36;
+  T x38 = T{0.0277777777777778}*u_2;
+  T x41 = x16 + u_2;
+  T x42 = T{6.00000000000000}*u_2;
+  T x43 = x18 + x42;
+  T x44 = rho*(x26 + x43 + T{9.00000000000000}*(x41*x41));
+  T x45 = -u_2;
+  T x46 = x45 + u_0;
+  T x47 = -x42;
+  T x48 = x18 + x47;
+  T x49 = rho*(x31 + x48 + T{9.00000000000000}*(x46*x46));
+  T x50 = u_0 + u_2;
+  T x51 = T{9.00000000000000}*(x50*x50);
+  T x52 = rho*(x26 + x48 + x51) + rho*(x31 + x43 + x51) - T{72.0000000000000}*f_curr[15] - T{72.0000000000000}*f_curr[3];
+  T x53 = x39 + x40 - x44 - x49 + x52;
+  T x54 = x13*x53;
+  T x57 = x25 + T{6.00000000000000}*x9;
+  T x58 = -x14 - x15 + x27 + x32 + x35;
+  T x59 = -x39 - x40 + x44 + x49 + x52;
+  T x60 = x56*(x21 + x57) + x56*(x22 + x57) + x58 + x59 - T{72.0000000000000}*f_curr[10] - T{72.0000000000000}*f_curr[8];
+  T x63 = x28 + u_2;
+  T x64 = x25 + x30;
+  T x65 = rho*(x42 + x64 + T{9.00000000000000}*(x63*x63));
+  T x66 = x45 + u_1;
+  T x67 = x20 + x25;
+  T x68 = rho*(x47 + x67 + T{9.00000000000000}*(x66*x66));
+  T x69 = u_1 + u_2;
+  T x70 = T{9.00000000000000}*(x69*x69);
+  T x71 = rho*(x42 + x67 + x70) + rho*(x47 + x64 + x70) - T{72.0000000000000}*f_curr[0] - T{72.0000000000000}*f_curr[18];
+  T x72 = x61 + x62 - x65 - x68 + x71;
+  T x73 = T{6.00000000000000}*x10 + x24;
+  T x74 = -x61 - x62 + x65 + x68 + x71;
+  T x75 = x56*(x20 + x73) + x56*(x30 + x73) + x58 + x74 - T{72.0000000000000}*f_curr[12] - T{72.0000000000000}*f_curr[6];
+  T x76 = T{6.00000000000000}*x11 + x23 + T{2.00000000000000};
+  T x77 = x56*(x43 + x76) + x56*(x48 + x76) + x59 + x74 - T{72.0000000000000}*f_curr[16] - T{72.0000000000000}*f_curr[2];
+  T x78 = ((x36*u_0 + x72*u_2 + x75*u_1)*u_1 + (x36*u_1 + x53*u_2 + x60*u_0)*u_0 + (x53*u_0 + x72*u_1 + x77*u_2)*u_2)/x12;
+  T x79 = x13*x72;
+  T n_0 = -x13*x55*x60 - x37*x8 - x38*x54 + x55*x78;
+  T n_1 = -x13*x75*x8 - x37*x55 - x38*x79 + x78*x8;
+  T n_2 = -x13*x38*x77 + x38*x78 - x54*x55 - x79*x8;
+
+  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;
+
+  float3 n = normalize(make_float3(n_0, n_1, n_2));
+  cell_shear_normal[3*gid+0] = n.x;
+  cell_shear_normal[3*gid+1] = n.y;
+  cell_shear_normal[3*gid+2] = n.z;
+}
+
+};
+
+struct CollectShearLayerVisibilityF {
+
+using call_tag = tag::post_process_by_spatial_cell_mask;
+
+template <typename T, typename S>
+__device__ static void apply(
+    descriptor::D3Q19
+  , std::size_t gid
+  , std::size_t iX
+  , std::size_t iY
+  , std::size_t iZ
+  , T* shear_normal
+  , float3 view_direction
+  , cudaSurfaceObject_t surface
+) {
+  float3 n = make_float3(shear_normal[3*gid+0], shear_normal[3*gid+1], shear_normal[3*gid+2]);
+  float visibility = dot(n, view_direction);
+  surf3Dwrite(visibility, surface, iX*sizeof(float), iY, iZ);
+}
+
+};
diff --git a/tangle/LLBM/kernel/collect_streamlines.h b/tangle/LLBM/kernel/collect_streamlines.h
new file mode 100644
index 0000000..c787713
--- /dev/null
+++ b/tangle/LLBM/kernel/collect_streamlines.h
@@ -0,0 +1,48 @@
+#pragma once
+
+__device__ void writeDot(cudaSurfaceObject_t texture, int x, int y, uchar4 pixel) {
+  surf2Dwrite(pixel, texture, (x-1)*sizeof(uchar4), (y  ));
+  surf2Dwrite(pixel, texture, (x-1)*sizeof(uchar4), (y+1));
+  surf2Dwrite(pixel, texture, (x  )*sizeof(uchar4), (y  ));
+  surf2Dwrite(pixel, texture, (x  )*sizeof(uchar4), (y+1));
+} 
+
+__device__ float2 readVelocity(cudaSurfaceObject_t u_x, cudaSurfaceObject_t u_y, float2 p) {
+  return make_float2(tex2D<float>(u_x, p.x, p.y), tex2D<float>(u_y, p.x, p.y));
+}
+
+template <typename DESCRIPTOR, typename T>
+__global__ void renderStreamlinesToTexture(
+    descriptor::Cuboid<DESCRIPTOR> cuboid
+  , cudaSurfaceObject_t u_x
+  , cudaSurfaceObject_t u_y
+  , T* origins
+  , unsigned nSteps
+  , float charU
+  , cudaSurfaceObject_t texture
+)  {
+  const int iOrigin = threadIdx.x + blockIdx.x * blockDim.x;
+
+  float h = 1;
+  float2 p = make_float2(origins[2*iOrigin+0], origins[2*iOrigin+1]);
+
+  for (unsigned iStep = 0; iStep < nSteps; ++iStep) {
+    float2 u = make_float2(tex2D<float>(u_x, p.x, p.y), tex2D<float>(u_y, p.x, p.y));
+
+    float2 k1 = readVelocity(u_x, u_y, p) / charU;
+    float2 k2 = readVelocity(u_x, u_y, p + 0.5f*h * k1) / charU;
+    float2 k3 = readVelocity(u_x, u_y, p + 0.5f*h * k2) / charU;
+    float2 k4 = readVelocity(u_x, u_y, p +      h * k3) / charU;
+
+    p += h * (1.f/6.f*k1 + 1.f/3.f*k2 + 1.f/3.f*k3 + 1.f/6.f*k4);
+
+    const int screenX = std::nearbyint(p.x);
+    const int screenY = cuboid.nY-1 - std::nearbyint(p.y);
+
+    if (screenX < 1 || screenY < 1 || screenX > cuboid.nX-2 || screenY > cuboid.nY-2) {
+      break;
+    }
+
+    writeDot(texture, screenX, screenY, {255,255,255,255});
+  }
+}
diff --git a/tangle/LLBM/kernel/collect_velocity_norm.h b/tangle/LLBM/kernel/collect_velocity_norm.h
new file mode 100644
index 0000000..ea099a8
--- /dev/null
+++ b/tangle/LLBM/kernel/collect_velocity_norm.h
@@ -0,0 +1,61 @@
+#pragma once
+#include <LLBM/call_tag.h>
+
+struct CollectVelocityNormF {
+
+using call_tag = tag::post_process_by_spatial_cell_mask;
+
+template <typename T, typename S>
+__device__ static void apply(
+    descriptor::D2Q9
+  , std::size_t gid
+  , std::size_t iX
+  , std::size_t iY
+  , std::size_t iZ
+  , T* u
+  , cudaSurfaceObject_t surface
+) {
+  float norm = length(make_float2(u[2*gid+0], u[2*gid+1]));
+  surf2Dwrite(norm, surface, iX*sizeof(float), iY);
+}
+
+template <typename T, typename S>
+__device__ static void apply(
+    descriptor::D3Q19
+  , std::size_t gid
+  , std::size_t iX
+  , std::size_t iY
+  , std::size_t iZ
+  , T* u
+  , cudaSurfaceObject_t surface
+  , T* u_norm = nullptr
+) {
+  float norm = length(make_float3(u[3*gid+0], u[3*gid+1], u[3*gid+2]));
+  surf3Dwrite(norm, surface, iX*sizeof(float), iY, iZ);
+  if (u_norm != nullptr) {
+    u_norm[gid] = norm;
+  }
+}
+
+};
+
+template <typename SLICE, typename SAMPLE, typename PALETTE>
+__global__ void renderSliceViewToTexture(std::size_t width, std::size_t height, SLICE slice, SAMPLE sample, PALETTE palette, cudaSurfaceObject_t texture) {
+  const int screenX = threadIdx.x + blockIdx.x * blockDim.x;
+  const int screenY = threadIdx.y + blockIdx.y * blockDim.y;
+
+  if (screenX > width-1 || screenY > height-1) {
+    return;
+  }
+  
+  const std::size_t gid = slice(screenX,screenY);
+  float3 color = palette(sample(gid));
+
+  uchar4 pixel {
+    static_cast<unsigned char>(color.x * 255),
+    static_cast<unsigned char>(color.y * 255),
+    static_cast<unsigned char>(color.z * 255),
+    255
+  };
+  surf2Dwrite(pixel, texture, screenX*sizeof(uchar4), screenY);
+}
diff --git a/tangle/LLBM/kernel/collide.h b/tangle/LLBM/kernel/collide.h
new file mode 100644
index 0000000..7ca02bf
--- /dev/null
+++ b/tangle/LLBM/kernel/collide.h
@@ -0,0 +1,131 @@
+#pragma once
+#include <LLBM/call_tag.h>
+
+struct BgkCollideO {
+
+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 tau) {
+  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]);
+  T x0 = T{1} / (tau);
+  T x1 = T{0.0138888888888889}*x0;
+  T x2 = T{6.00000000000000}*u_1;
+  T x3 = T{6.00000000000000}*u_0;
+  T x4 = u_0 + u_1;
+  T x5 = T{9.00000000000000}*(x4*x4);
+  T x6 = u_1*u_1;
+  T x7 = T{3.00000000000000}*x6;
+  T x8 = u_0*u_0;
+  T x9 = T{3.00000000000000}*x8;
+  T x10 = x9 + T{-2.00000000000000};
+  T x11 = x10 + x7;
+  T x12 = T{0.0555555555555556}*x0;
+  T x13 = -x3;
+  T x14 = T{2.00000000000000} - x7;
+  T x15 = x14 + T{6.00000000000000}*x8;
+  T x16 = -x9;
+  T x17 = x16 + x2;
+  T x18 = u_0 - u_1;
+  T x19 = x14 + T{9.00000000000000}*(x18*x18);
+  T x20 = T{6.00000000000000}*x6;
+  f_next[0] = -x1*(rho*(x11 + x2 + x3 - x5) + T{72.0000000000000}*f_curr[0]) + f_curr[0];
+  f_next[1] = x12*(rho*(x13 + x15) - T{18.0000000000000}*f_curr[1]) + f_curr[1];
+  f_next[2] = x1*(rho*(x13 + x17 + x19) - T{72.0000000000000}*f_curr[2]) + f_curr[2];
+  f_next[3] = -x12*(rho*(x10 + x2 - x20) + T{18.0000000000000}*f_curr[3]) + f_curr[3];
+  f_next[4] = -T{0.111111111111111}*x0*(T{2.00000000000000}*rho*x11 + T{9.00000000000000}*f_curr[4]) + f_curr[4];
+  f_next[5] = x12*(rho*(x17 + x20 + T{2.00000000000000}) - T{18.0000000000000}*f_curr[5]) + f_curr[5];
+  f_next[6] = x1*(rho*(x16 + x19 - x2 + x3) - T{72.0000000000000}*f_curr[6]) + f_curr[6];
+  f_next[7] = x12*(rho*(x15 + x3) - T{18.0000000000000}*f_curr[7]) + f_curr[7];
+  f_next[8] = x1*(rho*(x14 + x17 + x3 + x5) - T{72.0000000000000}*f_curr[8]) + f_curr[8];
+  
+}
+
+template <typename T, typename S>
+__device__ static void apply(descriptor::D3Q19, S f_curr[19], S f_next[19], std::size_t gid, T tau) {
+  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]);
+  T x0 = T{1} / (tau);
+  T x1 = T{0.0138888888888889}*x0;
+  T x2 = u_1 + u_2;
+  T x3 = T{9.00000000000000}*(x2*x2);
+  T x4 = T{6.00000000000000}*u_2;
+  T x5 = u_1*u_1;
+  T x6 = T{3.00000000000000}*x5;
+  T x7 = -x6;
+  T x8 = u_0*u_0;
+  T x9 = T{3.00000000000000}*x8;
+  T x10 = T{2.00000000000000} - x9;
+  T x11 = x10 + x7;
+  T x12 = x11 + x4;
+  T x13 = u_2*u_2;
+  T x14 = T{3.00000000000000}*x13;
+  T x15 = -x14;
+  T x16 = T{6.00000000000000}*u_1;
+  T x17 = x15 + x16;
+  T x18 = T{6.00000000000000}*u_0;
+  T x19 = -u_0;
+  T x20 = x19 + u_2;
+  T x21 = x14 + x6 + T{-2.00000000000000};
+  T x22 = x21 + x9;
+  T x23 = x22 - x4;
+  T x24 = T{0.0277777777777778}*x0;
+  T x25 = T{6.00000000000000}*x13;
+  T x26 = u_0 + u_2;
+  T x27 = T{9.00000000000000}*(x26*x26);
+  T x28 = x15 + x18;
+  T x29 = -u_1;
+  T x30 = x29 + u_2;
+  T x31 = x19 + u_1;
+  T x32 = -x16;
+  T x33 = x18 + x22;
+  T x34 = T{6.00000000000000}*x5;
+  T x35 = u_0 + u_1;
+  T x36 = T{9.00000000000000}*(x35*x35);
+  T x37 = T{6.00000000000000}*x8;
+  T x38 = x9 + T{-2.00000000000000};
+  T x39 = x29 + u_0;
+  T x40 = -x18 + x22;
+  T x41 = -u_2;
+  T x42 = x41 + u_1;
+  T x43 = x22 + x4;
+  T 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/tangle/LLBM/kernel/equilibrium_density_wall.h b/tangle/LLBM/kernel/equilibrium_density_wall.h
new file mode 100644
index 0000000..4f757f9
--- /dev/null
+++ b/tangle/LLBM/kernel/equilibrium_density_wall.h
@@ -0,0 +1,226 @@
+#pragma once
+#include <LLBM/call_tag.h>
+#include <LLBM/wall.h>
+#include <LLBM/descriptor.h>
+
+struct EquilibriumDensityWallO {
+
+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 rho_w, WallNormal<1,0>) {
+  T u = (rho_w - T{2.00000000000000}*f_curr[0] - T{2.00000000000000}*f_curr[1] - T{2.00000000000000}*f_curr[2] - f_curr[3] - f_curr[4] - f_curr[5])/rho_w;
+  T rho = rho_w;
+  T u_0 = u;
+  T u_1 = 0.;
+  T e0 = T{0.0277777777777778}*rho;
+  T e1 = T{3.00000000000000}*u_1;
+  T e2 = T{3.00000000000000}*u_0;
+  T e3 = u_0 + u_1;
+  T e4 = T{4.50000000000000}*(e3*e3);
+  T e5 = u_1*u_1;
+  T e6 = T{1.50000000000000}*e5;
+  T e7 = u_0*u_0;
+  T e8 = T{1.50000000000000}*e7;
+  T e9 = e8 + T{-1.00000000000000};
+  T e10 = e6 + e9;
+  T e11 = T{0.111111111111111}*rho;
+  T e12 = -e2;
+  T e13 = T{1.00000000000000} - e6;
+  T e14 = e13 + T{3.00000000000000}*e7;
+  T e15 = -e8;
+  T e16 = e1 + e15;
+  T e17 = u_0 - u_1;
+  T e18 = e13 + T{4.50000000000000}*(e17*e17);
+  T e19 = T{3.00000000000000}*e5;
+  f_next[0] = -e0*(e1 + e10 + e2 - e4);
+  f_next[1] = e11*(e12 + e14);
+  f_next[2] = e0*(e12 + e16 + e18);
+  f_next[3] = -e11*(e1 - e19 + e9);
+  f_next[4] = -T{0.444444444444444}*e10*rho;
+  f_next[5] = e11*(e16 + e19 + T{1.00000000000000});
+  f_next[6] = e0*(-e1 + e15 + e18 + e2);
+  f_next[7] = e11*(e14 + e2);
+  f_next[8] = e0*(e13 + e16 + e2 + e4);
+  
+}
+
+template <typename T, typename S>
+__device__ static void apply(descriptor::D2Q9, S f_curr[9], S f_next[9], std::size_t gid, T rho_w, WallNormal<-1,0>) {
+  T u = (-rho_w + f_curr[3] + f_curr[4] + f_curr[5] + T{2.00000000000000}*f_curr[6] + T{2.00000000000000}*f_curr[7] + T{2.00000000000000}*f_curr[8])/rho_w;
+  T rho = rho_w;
+  T u_0 = u;
+  T u_1 = 0.;
+  T e0 = T{0.0277777777777778}*rho;
+  T e1 = T{3.00000000000000}*u_1;
+  T e2 = T{3.00000000000000}*u_0;
+  T e3 = u_0 + u_1;
+  T e4 = T{4.50000000000000}*(e3*e3);
+  T e5 = u_1*u_1;
+  T e6 = T{1.50000000000000}*e5;
+  T e7 = u_0*u_0;
+  T e8 = T{1.50000000000000}*e7;
+  T e9 = e8 + T{-1.00000000000000};
+  T e10 = e6 + e9;
+  T e11 = T{0.111111111111111}*rho;
+  T e12 = -e2;
+  T e13 = T{1.00000000000000} - e6;
+  T e14 = e13 + T{3.00000000000000}*e7;
+  T e15 = -e8;
+  T e16 = e1 + e15;
+  T e17 = u_0 - u_1;
+  T e18 = e13 + T{4.50000000000000}*(e17*e17);
+  T e19 = T{3.00000000000000}*e5;
+  f_next[0] = -e0*(e1 + e10 + e2 - e4);
+  f_next[1] = e11*(e12 + e14);
+  f_next[2] = e0*(e12 + e16 + e18);
+  f_next[3] = -e11*(e1 - e19 + e9);
+  f_next[4] = -T{0.444444444444444}*e10*rho;
+  f_next[5] = e11*(e16 + e19 + T{1.00000000000000});
+  f_next[6] = e0*(-e1 + e15 + e18 + e2);
+  f_next[7] = e11*(e14 + e2);
+  f_next[8] = e0*(e13 + e16 + e2 + e4);
+  
+}
+
+template <typename T, typename S>
+__device__ static void apply(descriptor::D3Q19, S f_curr[19], S f_next[19], std::size_t gid, T rho_w, WallNormal<1,0,0>) {
+  T u = (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;
+  T rho = rho_w;
+  T u_0 = u;
+  T u_1 = 0.;
+  T u_2 = 0.;
+  T e0 = T{0.0277777777777778}*rho;
+  T e1 = u_1 + u_2;
+  T e2 = T{4.50000000000000}*(e1*e1);
+  T e3 = T{3.00000000000000}*u_2;
+  T e4 = u_1*u_1;
+  T e5 = T{1.50000000000000}*e4;
+  T e6 = -e5;
+  T e7 = u_0*u_0;
+  T e8 = T{1.50000000000000}*e7;
+  T e9 = T{1.00000000000000} - e8;
+  T e10 = e6 + e9;
+  T e11 = e10 + e3;
+  T e12 = T{3.00000000000000}*u_1;
+  T e13 = u_2*u_2;
+  T e14 = T{1.50000000000000}*e13;
+  T e15 = -e14;
+  T e16 = e12 + e15;
+  T e17 = T{3.00000000000000}*u_0;
+  T e18 = -u_2;
+  T e19 = e18 + u_0;
+  T e20 = -T{4.50000000000000}*e19*e19;
+  T e21 = e14 + e5 + T{-1.00000000000000};
+  T e22 = e21 + e8;
+  T e23 = e22 - e3;
+  T e24 = T{0.0555555555555556}*rho;
+  T e25 = T{3.00000000000000}*e13;
+  T e26 = u_0 + u_2;
+  T e27 = T{4.50000000000000}*(e26*e26);
+  T e28 = e15 + e17;
+  T e29 = e18 + u_1;
+  T e30 = -T{4.50000000000000}*e29*e29;
+  T e31 = -e12;
+  T e32 = u_0 - u_1;
+  T e33 = -T{4.50000000000000}*e32*e32;
+  T e34 = e17 + e22;
+  T e35 = T{3.00000000000000}*e4;
+  T e36 = u_0 + u_1;
+  T e37 = T{4.50000000000000}*(e36*e36);
+  T e38 = T{3.00000000000000}*e7;
+  T e39 = e8 + T{-1.00000000000000};
+  T e40 = -e17 + e22;
+  T 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 <typename T, typename S>
+__device__ static void apply(descriptor::D3Q19, S f_curr[19], S f_next[19], std::size_t gid, T rho_w, WallNormal<-1,0,0>) {
+  T u = (-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;
+  T rho = rho_w;
+  T u_0 = u;
+  T u_1 = 0.;
+  T u_2 = 0.;
+  T e0 = T{0.0277777777777778}*rho;
+  T e1 = u_1 + u_2;
+  T e2 = T{4.50000000000000}*(e1*e1);
+  T e3 = T{3.00000000000000}*u_2;
+  T e4 = u_1*u_1;
+  T e5 = T{1.50000000000000}*e4;
+  T e6 = -e5;
+  T e7 = u_0*u_0;
+  T e8 = T{1.50000000000000}*e7;
+  T e9 = T{1.00000000000000} - e8;
+  T e10 = e6 + e9;
+  T e11 = e10 + e3;
+  T e12 = T{3.00000000000000}*u_1;
+  T e13 = u_2*u_2;
+  T e14 = T{1.50000000000000}*e13;
+  T e15 = -e14;
+  T e16 = e12 + e15;
+  T e17 = T{3.00000000000000}*u_0;
+  T e18 = -u_2;
+  T e19 = e18 + u_0;
+  T e20 = -T{4.50000000000000}*e19*e19;
+  T e21 = e14 + e5 + T{-1.00000000000000};
+  T e22 = e21 + e8;
+  T e23 = e22 - e3;
+  T e24 = T{0.0555555555555556}*rho;
+  T e25 = T{3.00000000000000}*e13;
+  T e26 = u_0 + u_2;
+  T e27 = T{4.50000000000000}*(e26*e26);
+  T e28 = e15 + e17;
+  T e29 = e18 + u_1;
+  T e30 = -T{4.50000000000000}*e29*e29;
+  T e31 = -e12;
+  T e32 = u_0 - u_1;
+  T e33 = -T{4.50000000000000}*e32*e32;
+  T e34 = e17 + e22;
+  T e35 = T{3.00000000000000}*e4;
+  T e36 = u_0 + u_1;
+  T e37 = T{4.50000000000000}*(e36*e36);
+  T e38 = T{3.00000000000000}*e7;
+  T e39 = e8 + T{-1.00000000000000};
+  T e40 = -e17 + e22;
+  T 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/tangle/LLBM/kernel/equilibrium_velocity_wall.h b/tangle/LLBM/kernel/equilibrium_velocity_wall.h
new file mode 100644
index 0000000..23354dc
--- /dev/null
+++ b/tangle/LLBM/kernel/equilibrium_velocity_wall.h
@@ -0,0 +1,222 @@
+#pragma once
+#include <LLBM/call_tag.h>
+#include <LLBM/wall.h>
+#include <LLBM/descriptor.h>
+
+struct EquilibriumVelocityWallO {
+
+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_w, WallNormal<1,0>) {
+  T rho = -(T{2.00000000000000}*f_curr[0] + T{2.00000000000000}*f_curr[1] + T{2.00000000000000}*f_curr[2] + f_curr[3] + f_curr[4] + f_curr[5])/(u_w + T{-1.00000000000000});
+  T u_0 = u_w;
+  T u_1 = 0.;
+  T e0 = T{0.0277777777777778}*rho;
+  T e1 = T{3.00000000000000}*u_1;
+  T e2 = T{3.00000000000000}*u_0;
+  T e3 = u_0 + u_1;
+  T e4 = T{4.50000000000000}*(e3*e3);
+  T e5 = u_1*u_1;
+  T e6 = T{1.50000000000000}*e5;
+  T e7 = u_0*u_0;
+  T e8 = T{1.50000000000000}*e7;
+  T e9 = e8 + T{-1.00000000000000};
+  T e10 = e6 + e9;
+  T e11 = T{0.111111111111111}*rho;
+  T e12 = -e2;
+  T e13 = T{1.00000000000000} - e6;
+  T e14 = e13 + T{3.00000000000000}*e7;
+  T e15 = -e8;
+  T e16 = e1 + e15;
+  T e17 = u_0 - u_1;
+  T e18 = e13 + T{4.50000000000000}*(e17*e17);
+  T e19 = T{3.00000000000000}*e5;
+  f_next[0] = -e0*(e1 + e10 + e2 - e4);
+  f_next[1] = e11*(e12 + e14);
+  f_next[2] = e0*(e12 + e16 + e18);
+  f_next[3] = -e11*(e1 - e19 + e9);
+  f_next[4] = -T{0.444444444444444}*e10*rho;
+  f_next[5] = e11*(e16 + e19 + T{1.00000000000000});
+  f_next[6] = e0*(-e1 + e15 + e18 + e2);
+  f_next[7] = e11*(e14 + e2);
+  f_next[8] = e0*(e13 + e16 + e2 + e4);
+  
+}
+
+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_w, WallNormal<-1,0>) {
+  T rho = (f_curr[3] + f_curr[4] + f_curr[5] + T{2.00000000000000}*f_curr[6] + T{2.00000000000000}*f_curr[7] + T{2.00000000000000}*f_curr[8])/(u_w + T{1.00000000000000});
+  T u_0 = u_w;
+  T u_1 = 0;
+  T e0 = T{0.0277777777777778}*rho;
+  T e1 = T{3.00000000000000}*u_1;
+  T e2 = T{3.00000000000000}*u_0;
+  T e3 = u_0 + u_1;
+  T e4 = T{4.50000000000000}*(e3*e3);
+  T e5 = u_1*u_1;
+  T e6 = T{1.50000000000000}*e5;
+  T e7 = u_0*u_0;
+  T e8 = T{1.50000000000000}*e7;
+  T e9 = e8 + T{-1.00000000000000};
+  T e10 = e6 + e9;
+  T e11 = T{0.111111111111111}*rho;
+  T e12 = -e2;
+  T e13 = T{1.00000000000000} - e6;
+  T e14 = e13 + T{3.00000000000000}*e7;
+  T e15 = -e8;
+  T e16 = e1 + e15;
+  T e17 = u_0 - u_1;
+  T e18 = e13 + T{4.50000000000000}*(e17*e17);
+  T e19 = T{3.00000000000000}*e5;
+  f_next[0] = -e0*(e1 + e10 + e2 - e4);
+  f_next[1] = e11*(e12 + e14);
+  f_next[2] = e0*(e12 + e16 + e18);
+  f_next[3] = -e11*(e1 - e19 + e9);
+  f_next[4] = -T{0.444444444444444}*e10*rho;
+  f_next[5] = e11*(e16 + e19 + T{1.00000000000000});
+  f_next[6] = e0*(-e1 + e15 + e18 + e2);
+  f_next[7] = e11*(e14 + e2);
+  f_next[8] = e0*(e13 + e16 + e2 + e4);
+  
+}
+
+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_w, WallNormal<1,0,0>) {
+  T 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});
+  T u_0 = u_w;
+  T u_1 = 0;
+  T u_2 = 0;
+  T e0 = T{0.0277777777777778}*rho;
+  T e1 = u_1 + u_2;
+  T e2 = T{4.50000000000000}*(e1*e1);
+  T e3 = T{3.00000000000000}*u_2;
+  T e4 = u_1*u_1;
+  T e5 = T{1.50000000000000}*e4;
+  T e6 = -e5;
+  T e7 = u_0*u_0;
+  T e8 = T{1.50000000000000}*e7;
+  T e9 = T{1.00000000000000} - e8;
+  T e10 = e6 + e9;
+  T e11 = e10 + e3;
+  T e12 = T{3.00000000000000}*u_1;
+  T e13 = u_2*u_2;
+  T e14 = T{1.50000000000000}*e13;
+  T e15 = -e14;
+  T e16 = e12 + e15;
+  T e17 = T{3.00000000000000}*u_0;
+  T e18 = -u_2;
+  T e19 = e18 + u_0;
+  T e20 = -T{4.50000000000000}*e19*e19;
+  T e21 = e14 + e5 + T{-1.00000000000000};
+  T e22 = e21 + e8;
+  T e23 = e22 - e3;
+  T e24 = T{0.0555555555555556}*rho;
+  T e25 = T{3.00000000000000}*e13;
+  T e26 = u_0 + u_2;
+  T e27 = T{4.50000000000000}*(e26*e26);
+  T e28 = e15 + e17;
+  T e29 = e18 + u_1;
+  T e30 = -T{4.50000000000000}*e29*e29;
+  T e31 = -e12;
+  T e32 = u_0 - u_1;
+  T e33 = -T{4.50000000000000}*e32*e32;
+  T e34 = e17 + e22;
+  T e35 = T{3.00000000000000}*e4;
+  T e36 = u_0 + u_1;
+  T e37 = T{4.50000000000000}*(e36*e36);
+  T e38 = T{3.00000000000000}*e7;
+  T e39 = e8 + T{-1.00000000000000};
+  T e40 = -e17 + e22;
+  T 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 <typename T, typename S>
+__device__ static void apply(descriptor::D3Q19, S f_curr[19], S f_next[19], std::size_t gid, T u_w, WallNormal<-1,0,0>) {
+  T 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});
+  T u_0 = u_w;
+  T u_1 = 0;
+  T u_2 = 0;
+  T e0 = T{0.0277777777777778}*rho;
+  T e1 = u_1 + u_2;
+  T e2 = T{4.50000000000000}*(e1*e1);
+  T e3 = T{3.00000000000000}*u_2;
+  T e4 = u_1*u_1;
+  T e5 = T{1.50000000000000}*e4;
+  T e6 = -e5;
+  T e7 = u_0*u_0;
+  T e8 = T{1.50000000000000}*e7;
+  T e9 = T{1.00000000000000} - e8;
+  T e10 = e6 + e9;
+  T e11 = e10 + e3;
+  T e12 = T{3.00000000000000}*u_1;
+  T e13 = u_2*u_2;
+  T e14 = T{1.50000000000000}*e13;
+  T e15 = -e14;
+  T e16 = e12 + e15;
+  T e17 = T{3.00000000000000}*u_0;
+  T e18 = -u_2;
+  T e19 = e18 + u_0;
+  T e20 = -T{4.50000000000000}*e19*e19;
+  T e21 = e14 + e5 + T{-1.00000000000000};
+  T e22 = e21 + e8;
+  T e23 = e22 - e3;
+  T e24 = T{0.0555555555555556}*rho;
+  T e25 = T{3.00000000000000}*e13;
+  T e26 = u_0 + u_2;
+  T e27 = T{4.50000000000000}*(e26*e26);
+  T e28 = e15 + e17;
+  T e29 = e18 + u_1;
+  T e30 = -T{4.50000000000000}*e29*e29;
+  T e31 = -e12;
+  T e32 = u_0 - u_1;
+  T e33 = -T{4.50000000000000}*e32*e32;
+  T e34 = e17 + e22;
+  T e35 = T{3.00000000000000}*e4;
+  T e36 = u_0 + u_1;
+  T e37 = T{4.50000000000000}*(e36*e36);
+  T e38 = T{3.00000000000000}*e7;
+  T e39 = e8 + T{-1.00000000000000};
+  T e40 = -e17 + e22;
+  T 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/tangle/LLBM/kernel/executor.h b/tangle/LLBM/kernel/executor.h
new file mode 100644
index 0000000..942918d
--- /dev/null
+++ b/tangle/LLBM/kernel/executor.h
@@ -0,0 +1,171 @@
+#pragma once
+
+#include <LLBM/operator.h>
+
+namespace kernel {
+
+template <typename OPERATOR, typename DESCRIPTOR, typename T, typename S, typename... ARGS>
+__global__ void call_operator(
+    LatticeView<DESCRIPTOR,S> lattice
+  , std::size_t* cells
+  , std::size_t  cell_count
+  , ARGS... args
+) {
+  const std::size_t index = blockIdx.x * blockDim.x + threadIdx.x;
+  if (!(index < cell_count)) {
+      return;
+  }
+  const std::size_t gid = cells[index];
+
+  S f_curr[DESCRIPTOR::q];
+  S f_next[DESCRIPTOR::q];
+  S* preshifted_f[DESCRIPTOR::q];
+  for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+    preshifted_f[iPop] = lattice.pop(iPop, gid);
+    f_curr[iPop] = *preshifted_f[iPop];
+  }
+  OPERATOR::template apply<T,S>(DESCRIPTOR(), f_curr, f_next, gid, std::forward<ARGS>(args)...);
+  for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+    *preshifted_f[iPop] = f_next[iPop];
+  }
+}
+
+template <typename OPERATOR, typename DESCRIPTOR, typename T, typename S, typename... ARGS>
+__global__ void call_operator(
+    LatticeView<DESCRIPTOR,S> lattice
+  , bool* mask
+  , ARGS... args
+) {
+  const std::size_t gid = blockIdx.x * blockDim.x + threadIdx.x;
+  if (!(gid < lattice.cuboid.volume) || !mask[gid]) {
+      return;
+  }
+
+  S f_curr[DESCRIPTOR::q];
+  S f_next[DESCRIPTOR::q];
+  S* preshifted_f[DESCRIPTOR::q];
+  for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+    preshifted_f[iPop] = lattice.pop(iPop, gid);
+    f_curr[iPop] = *preshifted_f[iPop];
+  }
+  OPERATOR::template apply<T,S>(DESCRIPTOR(), f_curr, f_next, gid, std::forward<ARGS>(args)...);
+  for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+    *preshifted_f[iPop] = f_next[iPop];
+  }
+}
+
+template <typename FUNCTOR, typename DESCRIPTOR, typename T, typename S, typename... ARGS>
+__global__ void call_functor(
+    LatticeView<DESCRIPTOR,S> lattice
+  , bool* mask
+  , ARGS... args
+) {
+  const std::size_t gid = blockIdx.x * blockDim.x + threadIdx.x;
+  if (!(gid < lattice.cuboid.volume) || !mask[gid]) {
+      return;
+  }
+
+  S f_curr[DESCRIPTOR::q];
+  S* preshifted_f[DESCRIPTOR::q];
+  for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+    preshifted_f[iPop] = lattice.pop(iPop, gid);
+    f_curr[iPop] = *preshifted_f[iPop];
+  }
+  FUNCTOR::template apply<T,S>(DESCRIPTOR(), f_curr, gid, std::forward<ARGS>(args)...);
+}
+
+template <typename DESCRIPTOR, typename T, typename S, typename... OPERATOR>
+__global__ void call_operators(
+    LatticeView<DESCRIPTOR,S> lattice
+  , OPERATOR... ops
+) {
+  const std::size_t gid = blockIdx.x * blockDim.x + threadIdx.x;
+  if (!(gid < lattice.cuboid.volume)) {
+      return;
+  }
+
+  S f_curr[DESCRIPTOR::q];
+  S f_next[DESCRIPTOR::q];
+  S* preshifted_f[DESCRIPTOR::q];
+  for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+    preshifted_f[iPop] = lattice.pop(iPop, gid);
+    f_curr[iPop] = *preshifted_f[iPop];
+  }
+  (ops.template apply<DESCRIPTOR,T,S>(DESCRIPTOR(), f_curr, f_next, gid) || ... || false);
+  for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+    *preshifted_f[iPop] = f_next[iPop];
+  }
+}
+
+template <typename OPERATOR, typename DESCRIPTOR, typename T, typename S, typename... ARGS>
+__global__ void call_operator_using_list(
+    LatticeView<DESCRIPTOR,S> lattice
+  , std::size_t count
+  , ARGS... args
+) {
+  const std::size_t index = blockIdx.x * blockDim.x + threadIdx.x;
+  if (!(index < count)) {
+      return;
+  }
+  OPERATOR::template apply<T,S>(lattice, index, count, std::forward<ARGS>(args)...);
+}
+
+template <typename OPERATOR, typename DESCRIPTOR, typename T, typename S, typename... ARGS>
+__global__ void call_operator_using_list(
+    DESCRIPTOR descriptor
+  , std::size_t count
+  , ARGS... args
+) {
+  const std::size_t index = blockIdx.x * blockDim.x + threadIdx.x;
+  if (!(index < count)) {
+      return;
+  }
+  OPERATOR::template apply<T,S>(descriptor, index, count, std::forward<ARGS>(args)...);
+}
+
+template <typename FUNCTOR, typename DESCRIPTOR, typename T, typename S, typename... ARGS>
+__global__ void call_spatial_functor(
+    LatticeView<DESCRIPTOR,S> lattice
+  , bool* mask
+  , ARGS... args
+) {
+  const std::size_t iX = blockIdx.x * blockDim.x + threadIdx.x;
+  const std::size_t iY = blockIdx.y * blockDim.y + threadIdx.y;
+  const std::size_t iZ = blockIdx.z * blockDim.z + threadIdx.z;
+  if (!(iX < lattice.cuboid.nX && iY < lattice.cuboid.nY && iZ < lattice.cuboid.nZ)) {
+      return;
+  }
+  const std::size_t gid = descriptor::gid(lattice.cuboid,iX,iY,iZ);
+  if (!mask[gid]) {
+      return;
+  }
+
+  S f_curr[DESCRIPTOR::q];
+  S* preshifted_f[DESCRIPTOR::q];
+  for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+    preshifted_f[iPop] = lattice.pop(iPop, gid);
+    f_curr[iPop] = *preshifted_f[iPop];
+  }
+  FUNCTOR::template apply<T,S>(DESCRIPTOR(), f_curr, lattice.cuboid, gid, iX, iY, iZ, std::forward<ARGS>(args)...);
+}
+
+template <typename OPERATOR, typename DESCRIPTOR, typename T, typename S, typename... ARGS>
+__global__ void call_spatial_operator(
+    descriptor::Cuboid<DESCRIPTOR> cuboid
+  , bool* mask
+  , ARGS... args
+) {
+  const std::size_t iX = blockIdx.x * blockDim.x + threadIdx.x;
+  const std::size_t iY = blockIdx.y * blockDim.y + threadIdx.y;
+  const std::size_t iZ = blockIdx.z * blockDim.z + threadIdx.z;
+  if (!(iX < cuboid.nX && iY < cuboid.nY && iZ < cuboid.nZ)) {
+      return;
+  }
+  const std::size_t gid = descriptor::gid(cuboid,iX,iY,iZ);
+  if (!mask[gid]) {
+      return;
+  }
+  OPERATOR::template apply<T,S>(DESCRIPTOR(), gid, iX, iY, iZ, std::forward<ARGS>(args)...);
+}
+
+}
diff --git a/tangle/LLBM/kernel/free_slip.h b/tangle/LLBM/kernel/free_slip.h
new file mode 100644
index 0000000..25d7344
--- /dev/null
+++ b/tangle/LLBM/kernel/free_slip.h
@@ -0,0 +1,82 @@
+#pragma once
+#include <LLBM/call_tag.h>
+#include <LLBM/wall.h>
+#include <LLBM/descriptor.h>
+
+struct BounceBackFreeSlipO {
+
+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, WallNormal<1,0>) {
+  f_next[0] = f_curr[6];
+  f_next[1] = f_curr[7];
+  f_next[2] = f_curr[8];
+  f_next[3] = f_curr[3];
+  f_next[4] = f_curr[4];
+  f_next[5] = f_curr[5];
+  f_next[6] = f_curr[0];
+  f_next[7] = f_curr[1];
+  f_next[8] = f_curr[2];
+}
+
+template <typename T, typename S>
+__device__ static void apply(descriptor::D2Q9, S f_curr[9], S f_next[9], std::size_t gid, WallNormal<0,1>) {
+  f_next[0] = f_curr[2];
+  f_next[1] = f_curr[1];
+  f_next[2] = f_curr[0];
+  f_next[3] = f_curr[5];
+  f_next[4] = f_curr[4];
+  f_next[5] = f_curr[3];
+  f_next[6] = f_curr[8];
+  f_next[7] = f_curr[7];
+  f_next[8] = f_curr[6];
+}
+
+template <typename T, typename S>
+__device__ static void apply(descriptor::D3Q19, S f_curr[19], S f_next[19], std::size_t gid, WallNormal<0,1,0>) {
+  f_next[0] = f_curr[4];
+  f_next[1] = f_curr[1];
+  f_next[2] = f_curr[2];
+  f_next[3] = f_curr[3];
+  f_next[4] = f_curr[0];
+  f_next[5] = f_curr[11];
+  f_next[6] = f_curr[12];
+  f_next[7] = f_curr[13];
+  f_next[8] = f_curr[8];
+  f_next[9] = f_curr[9];
+  f_next[10] = f_curr[10];
+  f_next[11] = f_curr[5];
+  f_next[12] = f_curr[6];
+  f_next[13] = f_curr[7];
+  f_next[14] = f_curr[18];
+  f_next[15] = f_curr[15];
+  f_next[16] = f_curr[16];
+  f_next[17] = f_curr[17];
+  f_next[18] = f_curr[14];
+}
+
+template <typename T, typename S>
+__device__ static void apply(descriptor::D3Q19, S f_curr[19], S f_next[19], std::size_t gid, WallNormal<0,0,1>) {
+  f_next[0] = f_curr[14];
+  f_next[1] = f_curr[15];
+  f_next[2] = f_curr[16];
+  f_next[3] = f_curr[17];
+  f_next[4] = f_curr[18];
+  f_next[5] = f_curr[5];
+  f_next[6] = f_curr[6];
+  f_next[7] = f_curr[7];
+  f_next[8] = f_curr[8];
+  f_next[9] = f_curr[9];
+  f_next[10] = f_curr[10];
+  f_next[11] = f_curr[11];
+  f_next[12] = f_curr[12];
+  f_next[13] = f_curr[13];
+  f_next[14] = f_curr[0];
+  f_next[15] = f_curr[1];
+  f_next[16] = f_curr[2];
+  f_next[17] = f_curr[3];
+  f_next[18] = f_curr[4];
+}
+
+};
diff --git a/tangle/LLBM/kernel/initialize.h b/tangle/LLBM/kernel/initialize.h
new file mode 100644
index 0000000..221b9ad
--- /dev/null
+++ b/tangle/LLBM/kernel/initialize.h
@@ -0,0 +1,44 @@
+#pragma once
+#include <LLBM/call_tag.h>
+
+struct InitializeO {
+
+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] = T{0.0277777777777778};
+  f_next[1] = T{0.111111111111111};
+  f_next[2] = T{0.0277777777777778};
+  f_next[3] = T{0.111111111111111};
+  f_next[4] = T{0.444444444444444};
+  f_next[5] = T{0.111111111111111};
+  f_next[6] = T{0.0277777777777778};
+  f_next[7] = T{0.111111111111111};
+  f_next[8] = T{0.0277777777777778};
+}
+
+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] = 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};
+}
+
+};
diff --git a/tangle/LLBM/kernel/propagate.h b/tangle/LLBM/kernel/propagate.h
new file mode 100644
index 0000000..08f78ea
--- /dev/null
+++ b/tangle/LLBM/kernel/propagate.h
@@ -0,0 +1,20 @@
+#pragma once
+
+template <typename DESCRIPTOR, typename S>
+class LatticeView;
+
+template <typename DESCRIPTOR, typename S>
+__global__ void propagate(LatticeView<DESCRIPTOR,S> lattice, S** base, std::size_t size) {
+
+  for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+    std::ptrdiff_t shift = -descriptor::offset<DESCRIPTOR>(lattice.cuboid, iPop);
+
+    lattice.population[iPop] += shift;
+
+    if (lattice.population[iPop] < base[iPop]) {
+      lattice.population[iPop] += size;
+    } else if (lattice.population[iPop] + size > base[iPop] + 2*size) {
+      lattice.population[iPop] -= size;
+    }
+  }
+}
diff --git a/tangle/LLBM/kernel/smagorinsky_collide.h b/tangle/LLBM/kernel/smagorinsky_collide.h
new file mode 100644
index 0000000..3489479
--- /dev/null
+++ b/tangle/LLBM/kernel/smagorinsky_collide.h
@@ -0,0 +1,183 @@
+#pragma once
+#include <LLBM/call_tag.h>
+
+struct SmagorinskyBgkCollideO {
+
+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 tau, T smagorinsky) {
+  T x0 = f_curr[1] + f_curr[2];
+  T x1 = f_curr[3] + f_curr[6];
+  T x2 = x0 + x1 + f_curr[0] + f_curr[4] + f_curr[5] + f_curr[7] + f_curr[8];
+  T x3 = f_curr[0] - f_curr[8];
+  T x4 = T{1} / (x2);
+  T x5 = T{72.0000000000000}*f_curr[2];
+  T x6 = T{72.0000000000000}*f_curr[6];
+  T rho = x2;
+  T x31 = T{4.00000000000000}*rho;
+  T x40 = T{2.00000000000000}*rho;
+  T u_0 = -x4*(x0 + x3 - f_curr[6] - f_curr[7]);
+  T x7 = T{6.00000000000000}*u_0;
+  T x8 = -x7;
+  T x15 = u_0*u_0;
+  T x16 = T{3.00000000000000}*x15;
+  T x17 = -x16;
+  T x27 = x16 + T{-2.00000000000000};
+  T u_1 = -x4*(x1 + x3 - f_curr[2] - f_curr[5]);
+  T x9 = u_0 - u_1;
+  T x10 = T{9.00000000000000}*(x9*x9);
+  T x11 = u_1*u_1;
+  T x12 = T{3.00000000000000}*x11;
+  T x13 = T{2.00000000000000} - x12;
+  T x14 = T{6.00000000000000}*u_1;
+  T x18 = x14 + x17;
+  T x19 = x13 + x18;
+  T x20 = x10 + x19 + x8;
+  T x21 = rho*x20;
+  T x22 = x10 + x13 - x14 + x17 + x7;
+  T x23 = rho*x22;
+  T x24 = u_0 + u_1;
+  T x25 = T{9.00000000000000}*(x24*x24);
+  T x26 = x19 + x25 + x7;
+  T x28 = x12 + x27;
+  T x29 = x14 - x25 + x28 + x7;
+  T x30 = rho*x26 - rho*x29 - T{72.0000000000000}*f_curr[0] - T{72.0000000000000}*f_curr[8];
+  T x32 = x13 + T{6.00000000000000}*x15;
+  T x33 = x32 + x8;
+  T x34 = x32 + x7;
+  T x35 = x21 + x23 + x30 - x5 - x6;
+  T x36 = T{6.00000000000000}*x11;
+  T x37 = x14 + x27 - x36;
+  T x38 = x18 + x36 + T{2.00000000000000};
+  T x39 = T{1} / (tau + sqrt(T{0.707106781186548}*(smagorinsky*smagorinsky)*sqrt((-x21 - x23 + x30 + x5 + x6)*(-x21 - x23 + x30 + x5 + x6) + T{0.500000000000000}*((x31*x33 + x31*x34 + x35 - 72*f_curr[1] - 72*f_curr[7])*(x31*x33 + x31*x34 + x35 - 72*f_curr[1] - 72*f_curr[7])) + T{0.500000000000000}*((-x31*x37 + x31*x38 + x35 - 72*f_curr[3] - 72*f_curr[5])*(-x31*x37 + x31*x38 + x35 - 72*f_curr[3] - 72*f_curr[5]))) + tau*tau));
+  f_next[0] = -T{0.0138888888888889}*x39*(x29*x40 + T{144.000000000000}*f_curr[0]) + f_curr[0];
+  f_next[1] = T{0.0555555555555556}*x39*(x33*x40 - T{36.0000000000000}*f_curr[1]) + f_curr[1];
+  f_next[2] = T{0.0138888888888889}*x39*(x20*x40 - T{144.000000000000}*f_curr[2]) + f_curr[2];
+  f_next[3] = -T{0.0555555555555556}*x39*(x37*x40 + T{36.0000000000000}*f_curr[3]) + f_curr[3];
+  f_next[4] = -T{0.111111111111111}*x39*(T{4.00000000000000}*rho*x28 + T{18.0000000000000}*f_curr[4]) + f_curr[4];
+  f_next[5] = T{0.0555555555555556}*x39*(x38*x40 - T{36.0000000000000}*f_curr[5]) + f_curr[5];
+  f_next[6] = T{0.0138888888888889}*x39*(x22*x40 - T{144.000000000000}*f_curr[6]) + f_curr[6];
+  f_next[7] = T{0.0555555555555556}*x39*(x34*x40 - T{36.0000000000000}*f_curr[7]) + f_curr[7];
+  f_next[8] = T{0.0138888888888889}*x39*(x26*x40 - T{144.000000000000}*f_curr[8]) + f_curr[8];
+  
+}
+
+template <typename T, typename S>
+__device__ static void apply(descriptor::D3Q19, S f_curr[19], S f_next[19], std::size_t gid, T tau, T smagorinsky) {
+  T x0 = f_curr[10] + f_curr[13] + f_curr[17];
+  T x1 = f_curr[14] + f_curr[5] + f_curr[6];
+  T x2 = f_curr[1] + f_curr[2] + f_curr[4];
+  T 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];
+  T x4 = -f_curr[11] + f_curr[7];
+  T x5 = -f_curr[15] + f_curr[3];
+  T x6 = T{1} / (x3);
+  T x7 = f_curr[0] - f_curr[18];
+  T x8 = T{72.0000000000000}*f_curr[5];
+  T x9 = T{72.0000000000000}*f_curr[13];
+  T x42 = T{72.0000000000000}*f_curr[1];
+  T x43 = T{72.0000000000000}*f_curr[17];
+  T x61 = T{72.0000000000000}*f_curr[4];
+  T x62 = T{72.0000000000000}*f_curr[14];
+  T rho = x3;
+  T x74 = T{2.00000000000000}*rho;
+  T x89 = T{2.00000000000000}*rho;
+  T u_0 = x6*(x0 + x4 + x5 - f_curr[1] - f_curr[5] - f_curr[8]);
+  T x12 = -u_0;
+  T x14 = T{6.00000000000000}*u_0;
+  T x17 = u_0*u_0;
+  T x18 = T{3.00000000000000}*x17;
+  T x26 = -x14;
+  T x35 = T{2.00000000000000} - x18;
+  T x75 = T{6.00000000000000}*x17;
+  T u_1 = x6*(x1 + x4 + x7 - f_curr[12] - f_curr[13] - f_curr[4]);
+  T x10 = T{6.00000000000000}*u_1;
+  T x11 = -x10;
+  T x13 = x12 + u_1;
+  T x19 = u_1*u_1;
+  T x20 = T{3.00000000000000}*x19;
+  T x21 = x18 + x20 + T{-2.00000000000000};
+  T x27 = -u_1;
+  T x28 = x27 + u_0;
+  T x32 = u_0 + u_1;
+  T x33 = T{9.00000000000000}*(x32*x32);
+  T x34 = -x20;
+  T x36 = x34 + x35;
+  T x81 = T{6.00000000000000}*x19;
+  T u_2 = x6*(x2 + x5 + x7 - f_curr[14] - f_curr[16] - f_curr[17]);
+  T x15 = u_2*u_2;
+  T x16 = T{3.00000000000000}*x15;
+  T x22 = x16 + x21;
+  T x23 = x14 + x22;
+  T x24 = x11 + x23 - T{9.00000000000000}*x13*x13;
+  T x25 = rho*x24;
+  T x29 = x10 + x22;
+  T x30 = x26 + x29 - T{9.00000000000000}*x28*x28;
+  T x31 = rho*x30;
+  T x37 = -x16;
+  T x38 = x10 + x37;
+  T x39 = x14 + x33 + x36 + x38;
+  T x40 = x14 + x29 - x33;
+  T x41 = rho*x39 - rho*x40 - T{72.0000000000000}*f_curr[11] - T{72.0000000000000}*f_curr[7];
+  T x44 = T{6.00000000000000}*u_2;
+  T x45 = -x44;
+  T x46 = x12 + u_2;
+  T x47 = x23 + x45 - T{9.00000000000000}*x46*x46;
+  T x48 = rho*x47;
+  T x49 = -u_2;
+  T x50 = x49 + u_0;
+  T x51 = x22 + x44;
+  T x52 = x26 + x51 - T{9.00000000000000}*x50*x50;
+  T x53 = rho*x52;
+  T x54 = u_0 + u_2;
+  T x55 = T{9.00000000000000}*(x54*x54);
+  T x56 = x36 + x44;
+  T x57 = x14 + x37;
+  T x58 = x55 + x56 + x57;
+  T x59 = x14 + x51 - x55;
+  T x60 = rho*x58 - rho*x59 - T{72.0000000000000}*f_curr[15] - T{72.0000000000000}*f_curr[3];
+  T x63 = x27 + u_2;
+  T x64 = x29 + x45 - T{9.00000000000000}*x63*x63;
+  T x65 = rho*x64;
+  T x66 = x49 + u_1;
+  T x67 = x11 + x51 - T{9.00000000000000}*x66*x66;
+  T x68 = rho*x67;
+  T x69 = u_1 + u_2;
+  T x70 = T{9.00000000000000}*(x69*x69);
+  T x71 = x38 + x56 + x70;
+  T x72 = x29 + x44 - x70;
+  T x73 = rho*x71 - rho*x72 - T{72.0000000000000}*f_curr[0] - T{72.0000000000000}*f_curr[18];
+  T x76 = x16 + T{-2.00000000000000};
+  T x77 = x14 + x20 - x75 + x76;
+  T x78 = x34 + x57 + x75 + T{2.00000000000000};
+  T x79 = -x42 - x43 - x48 - x53 + x60;
+  T x80 = -x25 - x31 + x41 - x8 - x9;
+  T x82 = x10 + x18 + x76 - x81;
+  T x83 = x35 + x38 + x81;
+  T x84 = -x61 - x62 - x65 - x68 + x73;
+  T x85 = T{6.00000000000000}*x15;
+  T x86 = x21 + x44 - x85;
+  T x87 = x56 + x85;
+  T x88 = T{1} / (tau + sqrt(T{0.707106781186548}*(smagorinsky*smagorinsky)*sqrt((x25 + x31 + x41 + x8 + x9)*(x25 + x31 + x41 + x8 + x9) + (x42 + x43 + x48 + x53 + x60)*(x42 + x43 + x48 + x53 + x60) + (x61 + x62 + x65 + x68 + x73)*(x61 + x62 + x65 + x68 + x73) + T{0.500000000000000}*((-x74*x77 + x74*x78 + x79 + x80 - 72*f_curr[10] - 72*f_curr[8])*(-x74*x77 + x74*x78 + x79 + x80 - 72*f_curr[10] - 72*f_curr[8])) + T{0.500000000000000}*((-x74*x82 + x74*x83 + x80 + x84 - 72*f_curr[12] - 72*f_curr[6])*(-x74*x82 + x74*x83 + x80 + x84 - 72*f_curr[12] - 72*f_curr[6])) + T{0.500000000000000}*((-x74*x86 + x74*x87 + x79 + x84 - 72*f_curr[16] - 72*f_curr[2])*(-x74*x86 + x74*x87 + x79 + x84 - 72*f_curr[16] - 72*f_curr[2]))) + tau*tau));
+  f_next[0] = T{0.0138888888888889}*x88*(x71*x89 - T{144.000000000000}*f_curr[0]) + f_curr[0];
+  f_next[1] = -T{0.0138888888888889}*x88*(x47*x89 + T{144.000000000000}*f_curr[1]) + f_curr[1];
+  f_next[2] = T{0.0277777777777778}*x88*(x87*x89 - T{72.0000000000000}*f_curr[2]) + f_curr[2];
+  f_next[3] = T{0.0138888888888889}*x88*(x58*x89 - T{144.000000000000}*f_curr[3]) + f_curr[3];
+  f_next[4] = -T{0.0138888888888889}*x88*(x64*x89 + T{144.000000000000}*f_curr[4]) + f_curr[4];
+  f_next[5] = -T{0.0138888888888889}*x88*(x24*x89 + T{144.000000000000}*f_curr[5]) + f_curr[5];
+  f_next[6] = T{0.0277777777777778}*x88*(x83*x89 - T{72.0000000000000}*f_curr[6]) + f_curr[6];
+  f_next[7] = T{0.0138888888888889}*x88*(x39*x89 - T{144.000000000000}*f_curr[7]) + f_curr[7];
+  f_next[8] = -T{0.0277777777777778}*x88*(x77*x89 + T{72.0000000000000}*f_curr[8]) + f_curr[8];
+  f_next[9] = -T{0.166666666666667}*x88*(x22*x89 + T{12.0000000000000}*f_curr[9]) + f_curr[9];
+  f_next[10] = T{0.0277777777777778}*x88*(x78*x89 - T{72.0000000000000}*f_curr[10]) + f_curr[10];
+  f_next[11] = -T{0.0138888888888889}*x88*(x40*x89 + T{144.000000000000}*f_curr[11]) + f_curr[11];
+  f_next[12] = -T{0.0277777777777778}*x88*(x82*x89 + T{72.0000000000000}*f_curr[12]) + f_curr[12];
+  f_next[13] = -T{0.0138888888888889}*x88*(x30*x89 + T{144.000000000000}*f_curr[13]) + f_curr[13];
+  f_next[14] = -T{0.0138888888888889}*x88*(x67*x89 + T{144.000000000000}*f_curr[14]) + f_curr[14];
+  f_next[15] = -T{0.0138888888888889}*x88*(x59*x89 + T{144.000000000000}*f_curr[15]) + f_curr[15];
+  f_next[16] = -T{0.0277777777777778}*x88*(x86*x89 + T{72.0000000000000}*f_curr[16]) + f_curr[16];
+  f_next[17] = -T{0.0138888888888889}*x88*(x52*x89 + T{144.000000000000}*f_curr[17]) + f_curr[17];
+  f_next[18] = -T{0.0138888888888889}*x88*(x72*x89 + T{144.000000000000}*f_curr[18]) + f_curr[18];
+}
+
+};
diff --git a/tangle/LLBM/lattice.h b/tangle/LLBM/lattice.h
new file mode 100644
index 0000000..8ba1d66
--- /dev/null
+++ b/tangle/LLBM/lattice.h
@@ -0,0 +1,131 @@
+#pragma once
+
+#include "memory.h"
+#include "call_tag.h"
+#include "operator.h"
+
+#include "propagate.h"
+#include "kernel/initialize.h"
+#include "kernel/executor.h"
+
+template <typename DESCRIPTOR, typename T, typename S=T>
+class Lattice {
+private:
+const descriptor::Cuboid<DESCRIPTOR> _cuboid;
+
+CyclicPopulationBuffer<DESCRIPTOR,S> _population;
+
+public:
+Lattice(descriptor::Cuboid<DESCRIPTOR> cuboid):
+  _cuboid(cuboid),
+  _population(cuboid) { }
+
+descriptor::Cuboid<DESCRIPTOR> cuboid() const {
+  return _cuboid;
+}
+
+void stream() {
+  _population.stream();
+}
+
+template <typename... OPERATOR>
+void apply(OPERATOR... ops) {
+  const auto block_size = 32;
+  const auto block_count = (_cuboid.volume + block_size - 1) / block_size;
+  kernel::call_operators<DESCRIPTOR,T,S,OPERATOR...><<<block_count,block_size>>>(
+    _population.view(), ops...
+  );
+}
+
+template <typename OPERATOR, typename... ARGS>
+void apply(ARGS&&... args) {
+  call_operator<OPERATOR>(typename OPERATOR::call_tag{}, std::forward<ARGS&&>(args)...);
+}
+
+template <typename OPERATOR, typename... ARGS>
+void call_operator(tag::call_by_cell_id, DeviceBuffer<std::size_t>& cells, ARGS... args) {
+  const auto block_size = 32;
+  const auto block_count = (cells.size() + block_size - 1) / block_size;
+  kernel::call_operator<OPERATOR,DESCRIPTOR,T,S,ARGS...><<<block_count,block_size>>>(
+    _population.view(), cells.device(), cells.size(), std::forward<ARGS>(args)...
+  );
+}
+
+template <typename OPERATOR, typename... ARGS>
+void call_operator(tag::call_by_cell_id, DeviceBuffer<bool>& mask, ARGS... args) {
+  const auto block_size = 32;
+  const auto block_count = (_cuboid.volume + block_size - 1) / block_size;
+  kernel::call_operator<OPERATOR,DESCRIPTOR,T,S,ARGS...><<<block_count,block_size>>>(
+    _population.view(), mask.device(), std::forward<ARGS>(args)...
+  );
+}
+
+template <typename OPERATOR, typename... ARGS>
+void call_operator(tag::call_by_list_index, std::size_t count, ARGS... args) {
+  const auto block_size = 32;
+  const auto block_count = (count + block_size - 1) / block_size;
+  kernel::call_operator_using_list<OPERATOR,DESCRIPTOR,T,S,ARGS...><<<block_count,block_size>>>(
+    _population.view(), count, std::forward<ARGS>(args)...
+  );
+}
+
+template <typename FUNCTOR, typename... ARGS>
+void inspect(ARGS&&... args) {
+  call_functor<FUNCTOR>(typename FUNCTOR::call_tag{}, std::forward<ARGS&&>(args)...);
+}
+
+template <typename FUNCTOR, typename... ARGS>
+void call_functor(tag::call_by_cell_id, DeviceBuffer<std::size_t>& cells, ARGS... args) {
+  const auto block_size = 32;
+  const auto block_count = (cells.size() + block_size - 1) / block_size;
+  kernel::call_functor<FUNCTOR,DESCRIPTOR,T,S,ARGS...><<<block_count,block_size>>>(
+    _population.view(), cells.device(), cells.size(), std::forward<ARGS>(args)...
+  );
+}
+
+template <typename FUNCTOR, typename... ARGS>
+void call_functor(tag::call_by_cell_id, DeviceBuffer<bool>& mask, ARGS... args) {
+  const auto block_size = 32;
+  const auto block_count = (_cuboid.volume + block_size - 1) / block_size;
+  kernel::call_functor<FUNCTOR,DESCRIPTOR,T,S,ARGS...><<<block_count,block_size>>>(
+    _population.view(), mask.device(), std::forward<ARGS>(args)...
+  );
+}
+
+template <typename FUNCTOR, typename... ARGS>
+void call_functor(tag::call_by_spatial_cell_mask, DeviceBuffer<bool>& mask, ARGS... args) {
+  const dim3 block(32,8,4);
+  const dim3 grid((_cuboid.nX + block.x - 1) / block.x,
+                  (_cuboid.nY + block.y - 1) / block.y,
+                  (_cuboid.nZ + block.z - 1) / block.z);
+  kernel::call_spatial_functor<FUNCTOR,DESCRIPTOR,T,S,ARGS...><<<grid,block>>>(
+    _population.view(), mask.device(), std::forward<ARGS>(args)...
+  );
+}
+
+template <typename OPERATOR, typename... ARGS>
+void helper(ARGS&&... args) {
+  tagged_helper<OPERATOR>(typename OPERATOR::call_tag{}, std::forward<ARGS&&>(args)...);
+}
+
+template <typename OPERATOR, typename... ARGS>
+void tagged_helper(tag::post_process_by_list_index, std::size_t count, ARGS... args) {
+  const auto block_size = 32;
+  const auto block_count = (count + block_size - 1) / block_size;
+  kernel::call_operator_using_list<OPERATOR,DESCRIPTOR,T,S,ARGS...><<<block_count,block_size>>>(
+    DESCRIPTOR(), count, std::forward<ARGS>(args)...
+  );
+}
+
+template <typename OPERATOR, typename... ARGS>
+void tagged_helper(tag::post_process_by_spatial_cell_mask, DeviceBuffer<bool>& mask, ARGS... args) {
+  const dim3 block(32,8,4);
+  const dim3 grid((_cuboid.nX + block.x - 1) / block.x,
+                  (_cuboid.nY + block.y - 1) / block.y,
+                  (_cuboid.nZ + block.z - 1) / block.z);
+  kernel::call_spatial_operator<OPERATOR,DESCRIPTOR,T,S,ARGS...><<<grid,block>>>(
+    _cuboid, mask.device(), std::forward<ARGS>(args)...
+  );
+}
+
+};
diff --git a/tangle/LLBM/materials.h b/tangle/LLBM/materials.h
new file mode 100644
index 0000000..d782d8d
--- /dev/null
+++ b/tangle/LLBM/materials.h
@@ -0,0 +1,108 @@
+#pragma once
+
+#include "memory.h"
+#include "sdf.h"
+
+template <typename DESCRIPTOR>
+class CellMaterials : public SharedVector<int> {
+private:
+  const descriptor::Cuboid<DESCRIPTOR> _cuboid;
+  int* const _materials;
+
+public:
+  CellMaterials(descriptor::Cuboid<DESCRIPTOR> cuboid):
+    SharedVector<int>(cuboid.volume),
+    _cuboid(cuboid),
+    _materials(this->host()) { }
+
+  template <typename F>
+  CellMaterials(descriptor::Cuboid<DESCRIPTOR> cuboid, F f):
+    CellMaterials(cuboid) {
+    set(f);
+  }
+
+  descriptor::Cuboid<DESCRIPTOR> cuboid() const {
+    return _cuboid;
+  };
+ 
+  int get(std::size_t iCell) const {
+    return _materials[iCell];
+  }
+  
+  void set(std::size_t iCell, int material) {
+    _materials[iCell] = material;
+  }
+  template <typename F>
+  void set(F f) {
+    for (std::size_t iCell=0; iCell < _cuboid.volume; ++iCell) {
+      set(iCell, f(gidInverse(_cuboid, iCell)));
+    }
+  }
+  
+  template <typename S>
+  void sdf(S distance, int material, float eps=1e-2) {
+    for (std::size_t iCell=0; iCell < _cuboid.volume; ++iCell) {
+      auto p = gidInverseSmooth(_cuboid, iCell);
+      if (distance(p) < eps) {
+        set(iCell, material);
+      }
+    }
+  }
+  
+  void clean(int material) {
+    for (std::size_t iCell=0; iCell < _cuboid.volume; ++iCell) {
+      if (get(iCell) == material) {
+        if (_cuboid.isInside(iCell)) {
+          bool surrounded = true;
+          for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+            int m = get(descriptor::neighbor<DESCRIPTOR>(_cuboid, iCell, iPop));
+            surrounded &= m == material || m == 0;
+          }
+          if (surrounded) {
+            set(iCell, 0);
+          }
+        }
+      }
+    }
+  }
+  DeviceBuffer<std::size_t> list_of_material(int material) {
+  	std::vector<std::size_t> cells;
+  	for (std::size_t iCell=0; iCell < _cuboid.volume; ++iCell) {
+  		if (_materials[iCell] == material) {
+  			cells.emplace_back(iCell);
+  		}
+  	}
+  	return DeviceBuffer<std::size_t>(cells);
+  }
+  DeviceBuffer<bool> mask_of_material(int material) {
+    std::unique_ptr<bool[]> mask(new bool[_cuboid.volume]{});
+  	for (std::size_t iCell=0; iCell < _cuboid.volume; ++iCell) {
+      mask[iCell] = (_materials[iCell] == material);
+  	}
+  	return DeviceBuffer<bool>(mask.get(), _cuboid.volume);
+  }
+  std::size_t get_link_count(int bulk, int solid) {
+    std::size_t count = 0;
+    for (pop_index_t iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+      for (std::size_t iCell=0; iCell < _cuboid.volume; ++iCell) {
+        std::size_t jCell = descriptor::neighbor<DESCRIPTOR>(_cuboid, iCell, iPop);
+        if (get(iCell) == bulk && get(jCell) == solid) {
+          count++;
+        }
+      }
+    }
+    return count;
+  }
+  
+  template <typename F>
+  void for_links(int bulk, int solid, F f) {
+    for (pop_index_t iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+      for (std::size_t iCell=0; iCell < _cuboid.volume; ++iCell) {
+        std::size_t jCell = descriptor::neighbor<DESCRIPTOR>(_cuboid, iCell, iPop);
+        if (get(iCell) == bulk && get(jCell) == solid) {
+          f(iCell, iPop);
+        }
+      }
+    }
+  }
+};
diff --git a/tangle/LLBM/memory.h b/tangle/LLBM/memory.h
new file mode 100644
index 0000000..97fec5c
--- /dev/null
+++ b/tangle/LLBM/memory.h
@@ -0,0 +1,134 @@
+#pragma once
+
+#include <memory>
+#include <vector>
+#include <cstring>
+
+template <typename T>
+class DeviceBuffer {
+protected:
+  const std::size_t _size;
+  T* _data;
+
+public:
+  DeviceBuffer(std::size_t size):
+    _size(size) {
+    cudaMalloc(&_data, _size*sizeof(T));
+    cudaMemset(_data, 0, _size*sizeof(T));
+  }
+  DeviceBuffer(const T* data, std::size_t size):
+    DeviceBuffer(size) {
+    cudaMemcpy(_data, data, size*sizeof(T), cudaMemcpyHostToDevice);
+  }
+  DeviceBuffer(const std::vector<T>& data):
+    DeviceBuffer(data.data(), data.size()) { }
+ 
+  ~DeviceBuffer() {
+    cudaFree(_data);
+  }
+
+  T* device() {
+    return _data;
+  }
+
+  std::size_t size() const {
+    return _size;
+  }
+};
+
+template <typename T>
+class SharedVector : public DeviceBuffer<T> {
+private:
+  std::unique_ptr<T[]> _host_data;
+ 
+public:
+  SharedVector(std::size_t size):
+    DeviceBuffer<T>(size),
+    _host_data(new T[size]{}) {
+    syncDeviceFromHost();
+  }
+
+  T* host() {
+    return _host_data.get();
+  }
+
+  T& operator[](unsigned i) {
+    return host()[i];
+  }
+
+  void syncHostFromDevice() {
+    cudaMemcpy(_host_data.get(), this->_data, this->_size*sizeof(T), cudaMemcpyDeviceToHost);
+  }
+
+  void syncDeviceFromHost() {
+    cudaMemcpy(this->_data, _host_data.get(), this->_size*sizeof(T), cudaMemcpyHostToDevice);
+  }
+
+};
+
+template <typename T>
+class DeviceTexture {
+protected:
+  cudaExtent _extent;
+  cudaArray_t _array;
+
+  cudaChannelFormatDesc _channel_desc;
+  cudaResourceDesc _res_desc;
+  cudaTextureDesc  _tex_desc;
+
+  cudaTextureObject_t _texture;
+  cudaSurfaceObject_t _surface;
+
+public:
+  DeviceTexture(std::size_t nX, std::size_t nY, std::size_t nZ=0):
+    _extent(make_cudaExtent(nX,nY,nZ)),
+    _channel_desc(cudaCreateChannelDesc<float>()) {
+    cudaMalloc3DArray(&_array, &_channel_desc, _extent);
+
+    std::memset(&_res_desc, 0, sizeof(_res_desc));
+    _res_desc.resType = cudaResourceTypeArray;
+    _res_desc.res.array.array = _array;
+
+    std::memset(&_tex_desc, 0, sizeof(_tex_desc));
+    _res_desc.resType = cudaResourceTypeArray;
+    _tex_desc.addressMode[0]   = cudaAddressModeClamp;
+    _tex_desc.addressMode[1]   = cudaAddressModeClamp;
+    _tex_desc.addressMode[2]   = cudaAddressModeClamp;
+    _tex_desc.filterMode       = cudaFilterModeLinear;
+    _tex_desc.normalizedCoords = 0;
+
+    cudaCreateTextureObject(&_texture, &_res_desc, &_tex_desc, NULL);
+    cudaCreateSurfaceObject(&_surface, &_res_desc);
+  }
+
+  DeviceTexture(descriptor::CuboidD<3> c):
+    DeviceTexture(c.nX, c.nY, c.nZ) { }
+  
+  ~DeviceTexture() {
+    cudaFreeArray(_array);
+  }
+
+  cudaTextureObject_t getTexture() const {
+    return _texture;
+  }
+
+  cudaSurfaceObject_t getSurface() const {
+    return _surface;
+  }
+
+};
+
+__device__ float3 colorFromTexture(cudaSurfaceObject_t colormap, float value) {
+  uchar4 color{};
+  value = clamp(value, 0.f, 1.f);
+  surf2Dread(&color, colormap, unsigned(value * 999)*sizeof(uchar4), 0);
+  return make_float3(color.x / 255.f,
+                     color.y / 255.f,
+                     color.z / 255.f);
+}
+
+__device__ float noiseFromTexture(cudaSurfaceObject_t noisemap, int x, int y) {
+  uchar4 color{};
+  surf2Dread(&color, noisemap, x*sizeof(uchar4), y);
+  return color.x / 255.f;
+}
diff --git a/tangle/LLBM/operator.h b/tangle/LLBM/operator.h
new file mode 100644
index 0000000..6fb80e6
--- /dev/null
+++ b/tangle/LLBM/operator.h
@@ -0,0 +1,27 @@
+#pragma once
+
+#include <tuple>
+
+template <typename OPERATOR, typename... ARGS>
+struct Operator {
+  bool* const mask;
+  const std::tuple<ARGS...> config;
+
+  Operator(OPERATOR, DeviceBuffer<bool>& m, ARGS... args):
+    mask(m.device()),
+    config(args...) { }
+
+  template <typename DESCRIPTOR, typename T, typename S>
+  __device__ bool apply(DESCRIPTOR d, S f_curr[DESCRIPTOR::q], S f_next[DESCRIPTOR::q], std::size_t gid) const {
+    if (mask[gid]) {
+      std::apply([](auto... args) { OPERATOR::template apply<T,S>(args...); },
+                 std::tuple_cat(std::make_tuple(d, f_curr, f_next, gid), config));
+      return true;
+    } else {
+      return false;
+    }
+  }
+};
+
+template <typename OPERATOR, typename... ARGS>
+Operator(OPERATOR, DeviceBuffer<bool>&, ARGS... args) -> Operator<OPERATOR,std::remove_reference_t<ARGS>...>;
diff --git a/tangle/LLBM/propagate.h b/tangle/LLBM/propagate.h
new file mode 100644
index 0000000..d63ccd8
--- /dev/null
+++ b/tangle/LLBM/propagate.h
@@ -0,0 +1,111 @@
+#pragma once
+
+#include "memory.h"
+#include "descriptor.h"
+#include "kernel/propagate.h"
+
+#include <cuda.h>
+
+template <typename DESCRIPTOR, typename S>
+struct LatticeView {
+  const descriptor::Cuboid<DESCRIPTOR> cuboid;
+  S** population;
+
+  __device__ __forceinline__
+  S* pop(pop_index_t iPop, std::size_t gid) const;
+};
+
+template <typename DESCRIPTOR, typename S>
+class CyclicPopulationBuffer {
+protected:
+  const descriptor::Cuboid<DESCRIPTOR> _cuboid;
+
+  const std::size_t _page_size;
+  const std::size_t _volume;
+
+  CUmemGenericAllocationHandle _handle[DESCRIPTOR::q];
+  CUmemAllocationProp _prop{};
+  CUmemAccessDesc _access{};
+  CUdeviceptr _ptr;
+
+  SharedVector<S*> _base;
+  SharedVector<S*> _population;
+
+  S* device() {
+    return reinterpret_cast<S*>(_ptr);
+  }
+
+public:
+  CyclicPopulationBuffer(descriptor::Cuboid<DESCRIPTOR> cuboid);
+
+  LatticeView<DESCRIPTOR,S> view() {
+    return LatticeView<DESCRIPTOR,S>{ _cuboid, _population.device() };
+  }
+
+  void stream();
+
+};
+
+std::size_t getDevicePageSize(int device_id=-1) {
+  if (device_id == -1) {
+    cudaGetDevice(&device_id);
+  }
+  std::size_t granularity = 0;
+  CUmemAllocationProp prop = {};
+  prop.type = CU_MEM_ALLOCATION_TYPE_PINNED;
+  prop.location.type = CU_MEM_LOCATION_TYPE_DEVICE;
+  prop.location.id = device_id;
+  cuMemGetAllocationGranularity(&granularity, &prop, CU_MEM_ALLOC_GRANULARITY_MINIMUM);
+  return granularity;
+}
+
+template <typename DESCRIPTOR, typename S>
+CyclicPopulationBuffer<DESCRIPTOR,S>::CyclicPopulationBuffer(
+  descriptor::Cuboid<DESCRIPTOR> cuboid):
+  _cuboid(cuboid),
+  _page_size{getDevicePageSize()},
+  _volume{((cuboid.volume * sizeof(S) - 1) / _page_size + 1) * _page_size},
+  _base(DESCRIPTOR::q),
+  _population(DESCRIPTOR::q)
+{
+
+  int device_id = -1;
+  cudaGetDevice(&device_id);
+
+  _prop.type = CU_MEM_ALLOCATION_TYPE_PINNED;
+  _prop.location.type = CU_MEM_LOCATION_TYPE_DEVICE;
+  _prop.location.id = device_id;
+  cuMemAddressReserve(&_ptr, 2 * _volume * DESCRIPTOR::q, 0, 0, 0);
+
+  for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+    // per-population handle until cuMemMap accepts non-zero offset
+    cuMemCreate(&_handle[iPop], _volume, &_prop, 0); 
+    cuMemMap(_ptr + iPop * 2 * _volume,           _volume, 0, _handle[iPop], 0);
+    cuMemMap(_ptr + iPop * 2 * _volume + _volume, _volume, 0, _handle[iPop], 0);
+  }
+
+  _access.location.type = CU_MEM_LOCATION_TYPE_DEVICE;
+  _access.location.id = 0;
+  _access.flags = CU_MEM_ACCESS_FLAGS_PROT_READWRITE;
+  cuMemSetAccess(_ptr, 2 * _volume * DESCRIPTOR::q, &_access, 1);
+  cuMemsetD8(_ptr, 0, 2 * _volume * DESCRIPTOR::q);
+
+  for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+    _base[iPop] = device() + iPop * 2 * (_volume / sizeof(S));
+    _population[iPop] = _base[iPop] + iPop * ((_volume / sizeof(S)) / DESCRIPTOR::q);
+  }
+
+  _base.syncDeviceFromHost();
+  _population.syncDeviceFromHost();
+}
+
+template <typename DESCRIPTOR, typename S>
+__device__ __forceinline__
+S* LatticeView<DESCRIPTOR,S>::pop(pop_index_t iPop, std::size_t gid) const {
+  return population[iPop] + gid;
+}
+
+template <typename DESCRIPTOR, typename S>
+void CyclicPopulationBuffer<DESCRIPTOR,S>::stream() {
+  propagate<DESCRIPTOR,S><<<1,1>>>(view(), _base.device(), _volume / sizeof(S));
+}
diff --git a/tangle/LLBM/sdf.h b/tangle/LLBM/sdf.h
new file mode 100644
index 0000000..656c109
--- /dev/null
+++ b/tangle/LLBM/sdf.h
@@ -0,0 +1,75 @@
+#pragma once
+#include <vector_types.h>
+#include <cuda-samples/Common/helper_math.h>
+
+template <typename SDF, typename V>
+__device__ __host__
+float approximateDistance(SDF sdf, V origin, V dir, float d0, float d1, float eps=1e-2, unsigned N=128) {
+  float distance = d0;
+  float delta = (d1-d0) / N;
+  for (unsigned i=0; i < N; ++i) {
+    float d = sdf(origin + distance*dir);
+    if (d < eps) {
+      return distance;
+    }
+    distance += d;
+    if (distance > d1) {
+      return d1;
+    }
+  }
+  return d1;
+}
+
+namespace sdf {
+
+template <typename V>
+__device__ __host__ float sphere(V p, float r) {
+  return length(p) - r;
+}
+
+__device__ __host__ float box(float3 p, float3 b) {
+  float3 q = fabs(p) - b;
+  return length(fmaxf(q,make_float3(0))) + fmin(fmax(q.x,fmax(q.y,q.z)),0);
+}
+
+__device__ __host__ float cylinder(float3 p, float r, float h) {
+	return fmax(length(make_float2(p.x,p.y)) - r, fabs(p.z) - 0.5*h);
+}
+
+__device__ __host__ float add(float a, float b) {
+  return fmin(a, b);
+}
+
+__device__ __host__ float intersect(float a, float b) {
+  return fmax(a, b);
+}
+
+__device__ __host__ float sub(float a, float b) {
+  return intersect(-a, b);
+}
+
+__device__ __host__ float sadd(float a, float b, float k) {
+  float h = clamp(0.5f + 0.5f*(b-a)/k, 0.0f, 1.0f);
+  return lerp(b, a, h) - k*h*(1.f-h);
+}
+
+__device__ __host__ float ssub(float a, float b, float k) {
+  float h = clamp(0.5f - 0.5f*(b+a)/k, 0.f, 1.f);
+  return lerp(b, -a, h) + k*h*(1.f-h);
+}
+
+__device__ __host__ float sintersect(float a, float b, float k) {
+  float h = clamp(0.5f - 0.5f*(b-a)/k, 0.f, 1.f);
+  return lerp(b, a, h) + k*h*(1.f-h);
+}
+
+}
+
+template <typename SDF>
+__device__ float3 sdf_normal(SDF sdf, float3 v, float eps=1e-4) {
+  return normalize(make_float3(
+    sdf(make_float3(v.x + eps, v.y, v.z)) - sdf(make_float3(v.x - eps, v.y, v.z)),
+    sdf(make_float3(v.x, v.y + eps, v.z)) - sdf(make_float3(v.x, v.y - eps, v.z)),
+    sdf(make_float3(v.x, v.y, v.z + eps)) - sdf(make_float3(v.x, v.y, v.z - eps))
+  ));
+}
diff --git a/tangle/LLBM/sdf_boundary.h b/tangle/LLBM/sdf_boundary.h
new file mode 100644
index 0000000..f91480e
--- /dev/null
+++ b/tangle/LLBM/sdf_boundary.h
@@ -0,0 +1,114 @@
+#pragma once
+#include <LLBM/memory.h>
+#include <LLBM/materials.h>
+#include <LLBM/kernel/bouzidi.h>
+#include <iostream>
+
+template <typename DESCRIPTOR, typename T, typename S, typename SDF>
+class SignedDistanceBoundary {
+private:
+const descriptor::Cuboid<DESCRIPTOR> _cuboid;
+const std::size_t _count;
+
+SharedVector<std::size_t> _boundary;
+SharedVector<std::size_t> _fluid;
+SharedVector<std::size_t> _solid;
+SharedVector<S> _distance;
+SharedVector<S> _correction;
+SharedVector<S> _factor;
+SharedVector<pop_index_t> _missing;
+
+void set(std::size_t index, std::size_t iCell, pop_index_t iPop, S dist) {
+  pop_index_t jPop = descriptor::opposite<DESCRIPTOR>(iPop);
+  const std::size_t jPopCell = descriptor::neighbor<DESCRIPTOR>(_cuboid, iCell, jPop);
+  const std::size_t iPopCell = descriptor::neighbor<DESCRIPTOR>(_cuboid, iCell, iPop);
+
+  _boundary[index] = iCell;
+  _solid[index] = jPopCell;
+  _distance[index] = dist;
+  _correction[index] = 0;
+  _missing[index] = iPop;
+
+  T q = dist / descriptor::velocity_length<DESCRIPTOR>(iPop);
+  if (q > 0.5) {
+    _fluid[index] = iCell;
+    _factor[index] = 1 / (2*q);
+  } else {
+    _fluid[index] = iPopCell;
+    _factor[index] = 2*q;
+  }
+}
+
+void syncDeviceFromHost() {
+  _boundary.syncDeviceFromHost();
+  _fluid.syncDeviceFromHost();
+  _solid.syncDeviceFromHost();
+  _distance.syncDeviceFromHost();
+  _correction.syncDeviceFromHost();
+  _factor.syncDeviceFromHost();
+  _missing.syncDeviceFromHost();
+}
+
+public:
+SignedDistanceBoundary(Lattice<DESCRIPTOR,T,S>&, CellMaterials<DESCRIPTOR>& materials, SDF geometry, int bulk, int solid):
+  _cuboid(materials.cuboid()),
+  _count(materials.get_link_count(bulk, solid)),
+  _boundary(_count),
+  _fluid(_count),
+  _solid(_count),
+  _distance(_count),
+  _correction(_count),
+  _factor(_count),
+  _missing(_count)
+{
+  std::size_t index = 0;
+  materials.for_links(bulk, solid, [&](std::size_t iCell, pop_index_t iPop) {
+    auto p         = gidInverseSmooth(_cuboid, iCell);
+    auto direction = normalize(descriptor::velocity<DESCRIPTOR>(iPop));
+    float length   = descriptor::velocity_length<DESCRIPTOR>(iPop);
+    float distance = approximateDistance(geometry, p, direction, 0, length);
+    if (distance == 0.f || distance > length) {
+      std::cout << "Bogus distance d=" << distance << " at cell " << iCell
+                << " in direction " << std::to_string(iPop) << std::endl;
+    }
+    set(index++, iCell, descriptor::opposite<DESCRIPTOR>(iPop), distance);
+  });
+  syncDeviceFromHost();
+}
+
+template <typename VELOCITY>
+void setVelocity(VELOCITY field) {
+  for (std::size_t index=0; index < _count; ++index) {
+    pop_index_t jPop = descriptor::opposite<DESCRIPTOR>(_missing[index]);
+    auto direction = normalize(descriptor::velocity<DESCRIPTOR>(jPop));
+    float length = descriptor::velocity_length<DESCRIPTOR>(jPop);
+    auto p = descriptor::gidInverseSmooth(_cuboid, _boundary[index]);
+    auto u_w = field(p + _distance[index] * direction);
+    _correction[index] = 2*3*descriptor::weight<DESCRIPTOR>(jPop)
+                       * dot(u_w, descriptor::velocity<DESCRIPTOR>(jPop));
+    if (_distance[index] / length > 0.5) {
+      _correction[index] *= _factor[index];
+    }
+  }
+  _correction.syncDeviceFromHost();
+}
+
+std::size_t getCount() const {
+  return _count;
+}
+
+BouzidiConfig<S> getConfig() {
+  return BouzidiConfig<S>{
+    _boundary.device(),
+    _solid.device(),
+    _fluid.device(),
+    _factor.device(),
+    _correction.device(),
+    _missing.device()
+  };
+}
+
+};
+
+template <typename DESCRIPTOR, typename T, typename S, typename SDF>
+SignedDistanceBoundary(Lattice<DESCRIPTOR,T,S>&, CellMaterials<DESCRIPTOR>&, SDF, int, int) -> SignedDistanceBoundary<DESCRIPTOR,T,S,SDF>;
diff --git a/tangle/LLBM/volumetric.h b/tangle/LLBM/volumetric.h
new file mode 100644
index 0000000..312d0e8
--- /dev/null
+++ b/tangle/LLBM/volumetric.h
@@ -0,0 +1,134 @@
+#include <cuda-samples/Common/helper_math.h>
+
+#include <LLBM/sdf.h>
+
+__device__ float2 getNormalizedScreenPos(float w, float h, float x, float y) {
+	return make_float2(
+		2.f * (.5f - x/w) * w/h,
+		2.f * (.5f - y/h)
+	);
+}
+
+__device__ float3 getEyeRayDir(float2 screen_pos, float3 eye_pos, float3 eye_target) {
+	const float3 forward = normalize(eye_target - eye_pos);
+	const float3 right   = normalize(cross(make_float3(0.f, 0.f, -1.f), forward));
+	const float3 up      = normalize(cross(forward, right));
+
+	return normalize(screen_pos.x*right + screen_pos.y*up + 4*forward);
+}
+
+__device__ bool aabb(float3 origin, float3 dir, float3 min, float3 max, float& tmin, float& tmax) {
+  float3 invD = make_float3(1./dir.x, 1./dir.y, 1./dir.z);
+  float3 t0s = (min - origin) * invD;
+  float3 t1s = (max - origin) * invD;
+  float3 tsmaller = fminf(t0s, t1s);
+  float3 tbigger  = fmaxf(t0s, t1s);
+  tmin = fmaxf(tmin, fmaxf(tsmaller.x, fmaxf(tsmaller.y, tsmaller.z)));
+  tmax = fminf(tmax, fminf(tbigger.x, fminf(tbigger.y, tbigger.z)));
+  return (tmin < tmax);
+}
+
+__device__ bool aabb(float3 origin, float3 dir, descriptor::CuboidD<3>& cuboid, float& tmin, float& tmax) {
+  return aabb(origin, dir, make_float3(0), make_float3(cuboid.nX,cuboid.nY,cuboid.nZ), tmin, tmax);
+}
+
+struct VolumetricRenderConfig {
+  descriptor::CuboidD<3> cuboid;
+
+  cudaSurfaceObject_t palette;
+  cudaSurfaceObject_t noise;
+
+  float delta = 1;
+  float transparency = 1;
+  float brightness = 1;
+  float3 background = make_float3(22.f / 255.f);
+
+  float3 eye_pos;
+  float3 eye_dir;
+
+  cudaSurfaceObject_t canvas;
+  uint2 canvas_size;
+
+  bool align_slices_to_view = true;
+  bool apply_noise = true;
+  bool apply_blur = true;
+
+  VolumetricRenderConfig(descriptor::CuboidD<3> c):
+    cuboid(c) { }
+};
+
+
+
+template <typename SDF, typename SAMPLER, typename ATTENUATOR>
+__global__ void raymarch(
+  VolumetricRenderConfig config,
+  SDF geometry,
+  SAMPLER sampler,
+  ATTENUATOR attenuator
+) {
+  unsigned int x = blockIdx.x*blockDim.x + threadIdx.x;
+  unsigned int y = blockIdx.y*blockDim.y + threadIdx.y;
+
+  if (x > config.canvas_size.x - 1 || y > config.canvas_size.y - 1) {
+    return;
+  }
+
+  const float2 screen_pos = getNormalizedScreenPos(config.canvas_size.x, config.canvas_size.y, x, y);
+  const float3 ray_dir = getEyeRayDir(screen_pos, config.eye_pos, config.eye_pos + config.eye_dir);
+
+  float3 r = make_float3(0);
+  float  a = 0;
+  
+  float tmin = 0;
+  float tmax = 4000;
+  
+  if (aabb(config.eye_pos, ray_dir, config.cuboid, tmin, tmax)) {
+    float volume_dist = tmax - tmin;
+    float3 geometry_pos = config.eye_pos + tmin*ray_dir;
+    float geometry_dist = approximateDistance(geometry, geometry_pos, ray_dir, 0, volume_dist);
+    geometry_pos += geometry_dist * ray_dir;
+  
+    float jitter = config.align_slices_to_view * (floor(fabs(dot(config.eye_dir, tmin*ray_dir)) / config.delta) * config.delta - tmin)
+                 + config.apply_noise          * config.delta * noiseFromTexture(config.noise, threadIdx.x, threadIdx.y);
+  
+    tmin          += jitter;
+    volume_dist   -= jitter;
+    geometry_dist -= jitter;
+  
+    if (volume_dist > config.delta) {
+      float3 sample_pos = config.eye_pos + tmin * ray_dir;
+      unsigned n_samples = floor(geometry_dist / config.delta);
+      for (unsigned i=0; i < n_samples; ++i) {
+        sample_pos += config.delta * ray_dir;
+        
+        float  sample_value = sampler(sample_pos);
+        float3 sample_color = config.brightness * colorFromTexture(config.palette, sample_value);
+        
+        float sample_attenuation = attenuator(sample_value) * config.transparency;
+        float attenuation = 1 - a;
+        
+        r += attenuation * sample_attenuation * sample_color;
+        a += attenuation * sample_attenuation;
+      }
+    }
+  
+    if (geometry_dist < volume_dist) {
+      float3 n = sdf_normal(geometry, geometry_pos);
+      r = lerp((0.3f + fabs(dot(n, ray_dir))) * make_float3(0.3f), r, a);
+    }
+  } else {
+    a = 0;
+  }
+  
+  if (a < 1) {
+    r += (1 - a) * config.background;
+  }
+
+  uchar4 pixel {
+    static_cast<unsigned char>(clamp(r.x, 0.0f, 1.0f) * 255),
+    static_cast<unsigned char>(clamp(r.y, 0.0f, 1.0f) * 255),
+    static_cast<unsigned char>(clamp(r.z, 0.0f, 1.0f) * 255),
+    255
+  };
+  surf2Dwrite(pixel, config.canvas, x*sizeof(uchar4), y);
+}
diff --git a/tangle/LLBM/wall.h b/tangle/LLBM/wall.h
new file mode 100644
index 0000000..311fdd4
--- /dev/null
+++ b/tangle/LLBM/wall.h
@@ -0,0 +1,4 @@
+#pragma once
+
+template <int N_0, int N_1, int N_2=0>
+struct WallNormal { };
diff --git a/tangle/asset/noise/blue_0.png b/tangle/asset/noise/blue_0.png
new file mode 100644
index 0000000..d1c0534
--- /dev/null
+++ b/tangle/asset/noise/blue_0.png
diff --git a/tangle/asset/noise/blue_1.png b/tangle/asset/noise/blue_1.png
new file mode 100644
index 0000000..c32ee28
--- /dev/null
+++ b/tangle/asset/noise/blue_1.png
diff --git a/tangle/asset/noise/blue_2.png b/tangle/asset/noise/blue_2.png
new file mode 100644
index 0000000..0db58de
--- /dev/null
+++ b/tangle/asset/noise/blue_2.png
diff --git a/tangle/asset/noise/blue_3.png b/tangle/asset/noise/blue_3.png
new file mode 100644
index 0000000..e0d06cf
--- /dev/null
+++ b/tangle/asset/noise/blue_3.png
diff --git a/tangle/asset/noise/blue_4.png b/tangle/asset/noise/blue_4.png
new file mode 100644
index 0000000..28d1b80
--- /dev/null
+++ b/tangle/asset/noise/blue_4.png
diff --git a/tangle/asset/palette/4wave_ROTB.png b/tangle/asset/palette/4wave_ROTB.png
new file mode 100644
index 0000000..06c4d8d
--- /dev/null
+++ b/tangle/asset/palette/4wave_ROTB.png
diff --git a/tangle/asset/palette/4wave_equal.png b/tangle/asset/palette/4wave_equal.png
new file mode 100644
index 0000000..95e1fa4
--- /dev/null
+++ b/tangle/asset/palette/4wave_equal.png
diff --git a/tangle/asset/palette/5wave_cool.png b/tangle/asset/palette/5wave_cool.png
new file mode 100644
index 0000000..517896c
--- /dev/null
+++ b/tangle/asset/palette/5wave_cool.png
diff --git a/tangle/asset/palette/autumn.png b/tangle/asset/palette/autumn.png
new file mode 100644
index 0000000..e726a78
--- /dev/null
+++ b/tangle/asset/palette/autumn.png
diff --git a/tangle/asset/palette/blue.png b/tangle/asset/palette/blue.png
new file mode 100644
index 0000000..f2f6f3e
--- /dev/null
+++ b/tangle/asset/palette/blue.png
diff --git a/tangle/asset/palette/blue_orange.png b/tangle/asset/palette/blue_orange.png
new file mode 100644
index 0000000..a0b4ac2
--- /dev/null
+++ b/tangle/asset/palette/blue_orange.png
diff --git a/tangle/asset/palette/green_brown.png b/tangle/asset/palette/green_brown.png
new file mode 100644
index 0000000..ce23378
--- /dev/null
+++ b/tangle/asset/palette/green_brown.png
diff --git a/tangle/asset/palette/orange.png b/tangle/asset/palette/orange.png
new file mode 100644
index 0000000..c32f13b
--- /dev/null
+++ b/tangle/asset/palette/orange.png
diff --git a/tangle/asset/shader/blur.frag b/tangle/asset/shader/blur.frag
new file mode 100644
index 0000000..ad57bde
--- /dev/null
+++ b/tangle/asset/shader/blur.frag
@@ -0,0 +1,20 @@
+#version 330
+
+uniform sampler2D texture;
+
+layout(location = 0) out vec4 color;
+layout(origin_upper_left, pixel_center_integer) in vec4 gl_FragCoord;
+
+float kernel[7] = float[]( 0.00443184841193801, 0.0539909665131881, 0.241970724519143, 0.398942280401433, 0.241970724519143, 0.0539909665131881, 0.00443184841193801 );
+
+void main() {
+  vec3 blurred = vec3(0.0);
+
+  for (int i=-3; i <= 3; ++i) {
+    for (int j=-3; j <= 3; ++j) {
+      blurred += kernel[3+j] * kernel[3+i] * texelFetch(texture, ivec2(gl_FragCoord.xy) + ivec2(i,j), 0).xyz;
+    }
+  }
+
+  color = vec4(blurred, 1.0);
+}
diff --git a/tangle/benchmark-ldc.cu b/tangle/benchmark-ldc.cu
new file mode 100644
index 0000000..2443afe
--- /dev/null
+++ b/tangle/benchmark-ldc.cu
@@ -0,0 +1,95 @@
+#include <LLBM/base.h>
+
+#include <LLBM/kernel/collide.h>
+#include <LLBM/kernel/bounce_back.h>
+#include <LLBM/kernel/bounce_back_moving_wall.h>
+
+#include "util/timer.h"
+
+#include <iostream>
+
+using DESCRIPTOR = descriptor::D3Q19;
+
+template <typename T>
+void simulate(descriptor::Cuboid<DESCRIPTOR> cuboid, std::size_t nStep) {
+  cudaSetDevice(0);
+
+  Lattice<DESCRIPTOR,T> lattice(cuboid);
+  
+  CellMaterials<DESCRIPTOR> materials(cuboid, [&cuboid](uint3 p) -> int {
+    if (p.x == 0 || p.x == cuboid.nX-1 || p.y == 0 || p.y == cuboid.nY-1 || p.z == 0) {
+      return 2; // boundary cell
+    } else if (p.z == cuboid.nZ-1) {
+      return 3; // lid cell
+    } else {
+      return 1; // bulk
+    }
+  });
+  
+  auto bulk_mask = materials.mask_of_material(1);
+  auto box_mask  = materials.mask_of_material(2);
+  auto lid_mask  = materials.mask_of_material(3);
+  
+  auto bulk_cells = materials.list_of_material(1);
+  auto box_cells  = materials.list_of_material(2);
+  auto lid_cells  = materials.list_of_material(3);
+  
+  lattice.template apply<InitializeO>(bulk_cells);
+  lattice.template apply<InitializeO>(box_cells);
+  lattice.template apply<InitializeO>(lid_cells);
+  
+  cudaDeviceSynchronize();
+
+  for (std::size_t iStep=0; iStep < 100; ++iStep) {
+    lattice.apply(Operator(BgkCollideO(), bulk_mask, 0.56),
+                  Operator(BounceBackO(), box_mask),
+                  Operator(BounceBackMovingWallO(), lid_mask, 0.05f, 0.f, 0.f));
+    lattice.stream();
+  }
+
+  cudaDeviceSynchronize();
+
+  auto start = timer::now();
+
+  for (std::size_t iStep=0; iStep < nStep; ++iStep) {
+    lattice.apply(Operator(BgkCollideO(), bulk_mask, 0.56),
+                  Operator(BounceBackO(), box_mask),
+                  Operator(BounceBackMovingWallO(), lid_mask, 0.05f, 0.f, 0.f));
+    lattice.stream();
+  }
+
+  cudaDeviceSynchronize();
+
+  auto mlups = timer::mlups(cuboid.volume, nStep, start);
+
+  std::cout << sizeof(T) << ", " << cuboid.nX << ", " << nStep << ", " << mlups << std::endl;
+}
+
+int main(int argc, char* argv[]) {
+  if (argc < 3 || argc > 4) {
+    std::cerr << "Invalid parameter count" << std::endl;
+    return -1;
+  }
+
+  const std::size_t n     = atoi(argv[1]);
+  const std::size_t steps = atoi(argv[2]);
+
+  unsigned precision = 4;
+  if (argc == 4) {
+    precision = atoi(argv[3]);
+  }
+
+  switch (precision) {
+  case 4:
+    simulate<float>({ n, n, n}, steps);
+    break;
+  case 8:
+    simulate<double>({ n, n, n}, steps);
+    break;
+  default:
+    std::cerr << "Invalid precision" << std::endl;
+    return -1;
+  }
+
+  return 0;
+}
diff --git a/tangle/channel-with-sphere.cu b/tangle/channel-with-sphere.cu
new file mode 100644
index 0000000..d2effb0
--- /dev/null
+++ b/tangle/channel-with-sphere.cu
@@ -0,0 +1,85 @@
+#include <LLBM/base.h>
+#include <LLBM/bulk.h>
+#include <LLBM/boundary.h>
+
+#include "util/render_window.h"
+#include "util/texture.h"
+#include "util/colormap.h"
+
+#include "util/volumetric_example.h"
+#include "sampler/velocity_norm.h"
+#include "sampler/curl_norm.h"
+#include "sampler/q_criterion.h"
+
+using T = float;
+using DESCRIPTOR = descriptor::D3Q19;
+
+int main() {
+cudaSetDevice(0);
+
+const descriptor::Cuboid<DESCRIPTOR> cuboid(300, 80, 80);
+Lattice<DESCRIPTOR,T> lattice(cuboid);
+
+CellMaterials<DESCRIPTOR> materials(cuboid, [&cuboid](uint3 p) -> int {
+  if (p.z == 0 || p.z == cuboid.nZ-1) {
+    return 2; // boundary cell
+  } else if (p.y == 0 || p.y == cuboid.nY-1) {
+    return 3; // boundary cell
+  } else if (p.x == 0) {
+    return 4; // inflow cell
+  } else if (p.x == cuboid.nX-1) {
+    return 5; // outflow cell
+  } else {
+    return 1; // bulk
+  }
+});
+
+for (std::size_t iX=0; iX < cuboid.nX; ++iX) {
+  materials.set(gid(cuboid, iX, 0,           0), 6);
+  materials.set(gid(cuboid, iX, cuboid.nY-1, 0), 6);
+  materials.set(gid(cuboid, iX, 0,           cuboid.nZ-1), 6);
+  materials.set(gid(cuboid, iX, cuboid.nY-1, cuboid.nZ-1), 6);
+}
+
+auto obstacle = [cuboid] __host__ __device__ (float3 p) -> float {
+                  float3 q = p - make_float3(cuboid.nX/6, cuboid.nY/2, cuboid.nZ/2);
+                  return sdf::sphere(q, cuboid.nY/T{5});
+                };
+materials.sdf(obstacle, 0);
+SignedDistanceBoundary bouzidi(lattice, materials, obstacle, 1, 0);
+
+auto bulk_mask    = materials.mask_of_material(1);
+auto wall_mask_z  = materials.mask_of_material(2);
+auto wall_mask_y  = materials.mask_of_material(3);
+auto inflow_mask  = materials.mask_of_material(4);
+auto outflow_mask = materials.mask_of_material(5);
+auto edge_mask    = materials.mask_of_material(6);
+
+lattice.apply(Operator(InitializeO(), bulk_mask),
+              Operator(InitializeO(), wall_mask_z),
+              Operator(InitializeO(), wall_mask_y),
+              Operator(InitializeO(), inflow_mask),
+              Operator(InitializeO(), outflow_mask),
+              Operator(InitializeO(), edge_mask));
+
+cudaDeviceSynchronize();
+
+VolumetricExample renderer(cuboid);
+renderer.add<QCriterionS>(lattice, bulk_mask, obstacle);
+renderer.add<CurlNormS>(lattice, bulk_mask, obstacle);
+renderer.add<VelocityNormS>(lattice, bulk_mask, obstacle);
+renderer.run([&](std::size_t iStep) {
+  const float tau = 0.51;
+  const float inflow = 0.08;
+  
+  lattice.apply(Operator(BgkCollideO(), bulk_mask, tau),
+                Operator(BounceBackFreeSlipO(), wall_mask_z, WallNormal<0,0,1>()),
+                Operator(BounceBackFreeSlipO(), wall_mask_y, WallNormal<0,1,0>()),
+                Operator(EquilibriumVelocityWallO(), inflow_mask, std::min(iStep*1e-4, 1.0)*inflow, WallNormal<1,0,0>()),
+                Operator(EquilibriumDensityWallO(), outflow_mask, 1, WallNormal<-1,0,0>()),
+                Operator(BounceBackO(), edge_mask));
+  lattice.apply<BouzidiO>(bouzidi.getCount(), bouzidi.getConfig());
+  
+  lattice.stream();
+});
+}
diff --git a/tangle/ldc-2d.cu b/tangle/ldc-2d.cu
new file mode 100644
index 0000000..8989374
--- /dev/null
+++ b/tangle/ldc-2d.cu
@@ -0,0 +1,82 @@
+#include <LLBM/base.h>
+#include <LLBM/bulk.h>
+#include <LLBM/boundary.h>
+
+#include "util/render_window.h"
+#include "util/texture.h"
+#include "util/colormap.h"
+
+#include <LLBM/kernel/collect_moments.h>
+#include <LLBM/kernel/collect_velocity_norm.h>
+
+using T = float;
+using DESCRIPTOR = descriptor::D2Q9;
+
+int main() {
+cudaSetDevice(0);
+
+const descriptor::Cuboid<DESCRIPTOR> cuboid(500, 500);
+Lattice<DESCRIPTOR,T> lattice(cuboid);
+
+CellMaterials<DESCRIPTOR> materials(cuboid, [&cuboid](uint2 p) -> int {
+  if (p.x == 0 || p.y == 0 || p.x == cuboid.nX-1) {
+    return 2; // boundary cell
+  } else if (p.y == cuboid.nY-1) {
+    return 3; // lid cell
+  } else {
+    return 1; // bulk
+  }
+});
+
+auto bulk_mask = materials.mask_of_material(1);
+auto wall_mask = materials.mask_of_material(2);
+auto lid_mask  = materials.mask_of_material(3);
+
+lattice.apply(Operator(InitializeO(), bulk_mask),
+              Operator(InitializeO(), wall_mask),
+              Operator(InitializeO(), lid_mask));
+cudaDeviceSynchronize();
+
+const float tau = 0.51;
+const float u_lid = 0.05;
+
+RenderWindow window("LDC");
+cudaSurfaceObject_t colormap;
+ColorPalette palette(colormap);
+auto slice = [cuboid] __device__ (int iX, int iY) -> std::size_t {
+               return descriptor::gid(cuboid,iX,cuboid.nY-1-iY);
+             };
+DeviceBuffer<T> moments_rho(cuboid.volume);
+DeviceBuffer<T> moments_u(2*cuboid.volume);
+T* u = moments_u.device();
+std::size_t iStep = 0;
+
+while (window.isOpen()) {
+  lattice.apply(Operator(BgkCollideO(), bulk_mask, tau),
+                Operator(BounceBackO(), wall_mask),
+                Operator(BounceBackMovingWallO(), lid_mask, std::min(iStep*1e-3, 1.0)*u_lid, 0.f));
+  lattice.stream();
+  if (iStep % 100 == 0) {
+    cudaDeviceSynchronize();
+    lattice.inspect<CollectMomentsF>(bulk_mask, moments_rho.device(), moments_u.device());
+    renderSliceViewToTexture<<<
+      dim3(cuboid.nX / 32 + 1, cuboid.nY / 32 + 1),
+      dim3(32,32)
+    >>>(cuboid.nX, cuboid.nY,
+        slice,
+        [u,u_lid] __device__ (std::size_t gid) -> float {
+          return length(make_float2(u[2*gid+0], u[2*gid+1])) / u_lid;
+        },
+        [colormap] __device__ (float x) -> float3 {
+          return colorFromTexture(colormap, clamp(x, 0.f, 1.f));
+        },
+        window.getRenderSurface());
+    window.draw([&]() {
+      ImGui::Begin("Render");
+      palette.interact();
+      ImGui::End();
+    }, [](sf::Event&) { });
+  }
+  ++iStep;
+}
+}
diff --git a/tangle/ldc-3d.cu b/tangle/ldc-3d.cu
new file mode 100644
index 0000000..ece1234
--- /dev/null
+++ b/tangle/ldc-3d.cu
@@ -0,0 +1,58 @@
+#include <LLBM/base.h>
+#include <LLBM/bulk.h>
+#include <LLBM/boundary.h>
+
+#include "util/render_window.h"
+#include "util/texture.h"
+#include "util/colormap.h"
+
+#include "util/volumetric_example.h"
+#include "sampler/velocity_norm.h"
+#include "sampler/curl_norm.h"
+#include "sampler/shear_layer.h"
+
+using T = float;
+using DESCRIPTOR = descriptor::D3Q19;
+
+int main() {
+cudaSetDevice(0);
+
+const descriptor::Cuboid<DESCRIPTOR> cuboid(100, 100, 100);
+Lattice<DESCRIPTOR,T> lattice(cuboid);
+
+CellMaterials<DESCRIPTOR> materials(cuboid, [&cuboid](uint3 p) -> int {
+  if (p.x == 0 || p.x == cuboid.nX-1 || p.y == 0 || p.y == cuboid.nY-1 || p.z == cuboid.nZ-1) {
+    return 2; // boundary cell
+  } else if (p.z == 0) {
+    return 3; // lid cell
+  } else {
+    return 1; // bulk
+  }
+});
+
+auto bulk_mask = materials.mask_of_material(1);
+auto wall_mask = materials.mask_of_material(2);
+auto lid_mask  = materials.mask_of_material(3);
+
+lattice.apply(Operator(InitializeO(), bulk_mask),
+              Operator(InitializeO(), wall_mask),
+              Operator(InitializeO(), lid_mask));
+
+cudaDeviceSynchronize();
+
+auto none = [] __device__ (float3) -> float { return 1; };
+VolumetricExample renderer(cuboid);
+renderer.add<CurlNormS>(lattice, bulk_mask, none);
+renderer.add<ShearLayerVisibilityS>(lattice, bulk_mask, none, make_float3(0,1,0));
+renderer.add<VelocityNormS>(lattice, bulk_mask, none);
+renderer.run([&](std::size_t iStep) {
+  const float tau = 0.56;
+  const float lid = 0.10;
+  
+  lattice.apply(Operator(BgkCollideO(), bulk_mask, tau),
+                Operator(BounceBackO(), wall_mask),
+                Operator(BounceBackMovingWallO(), lid_mask, std::min(iStep*1e-3, 1.0)*lid, 0.f, 0.f));
+  
+  lattice.stream();
+});
+}
diff --git a/tangle/magnus.cu b/tangle/magnus.cu
new file mode 100644
index 0000000..aa31ba1
--- /dev/null
+++ b/tangle/magnus.cu
@@ -0,0 +1,116 @@
+#include <LLBM/base.h>
+#include <LLBM/bulk.h>
+#include <LLBM/boundary.h>
+
+#include "util/render_window.h"
+#include "util/texture.h"
+#include "util/colormap.h"
+
+#include <LLBM/kernel/collect_moments.h>
+#include <LLBM/kernel/collect_velocity_norm.h>
+
+using T = float;
+using DESCRIPTOR = descriptor::D2Q9;
+
+int main() {
+cudaSetDevice(0);
+
+const descriptor::Cuboid<DESCRIPTOR> cuboid(1200, 500);
+Lattice<DESCRIPTOR,T> lattice(cuboid);
+
+const float tau = 0.54;
+const float u_inflow = 0.02;
+const float u_rotate = 0.08;
+
+CellMaterials<DESCRIPTOR> materials(cuboid, [&cuboid](uint2 p) -> int {
+  if (p.x == 0) {
+    return 3; // inflow
+  } else if (p.x == cuboid.nX-1) {
+    return 4; // outflow
+  } else if (p.y == 0 || p.y == cuboid.nY-1) {
+    return 2; // wall
+  } else {
+    return 1; // bulk
+  }
+});
+
+materials.set(gid(cuboid, 0,0), 2);
+materials.set(gid(cuboid, 0,cuboid.nY-1), 2);
+materials.set(gid(cuboid, cuboid.nX-1,0), 5);
+materials.set(gid(cuboid, cuboid.nX-1,cuboid.nY-1), 5);
+
+auto cylinder = [cuboid] __host__ __device__ (float2 p) -> float {
+                  float2 q = p - make_float2(cuboid.nX/6, 3*cuboid.nY/4);
+                  float2 r = p - make_float2(cuboid.nX/6, 1*cuboid.nY/4);
+                  return sdf::add(sdf::sphere(q, cuboid.nY/18),
+                                  sdf::sphere(r, cuboid.nY/18));
+                };
+
+materials.sdf(cylinder, 0);
+SignedDistanceBoundary bouzidi(lattice, materials, cylinder, 1, 0);
+
+bouzidi.setVelocity([cuboid,u_rotate](float2 p) -> float2 {
+  float2 q = p - make_float2(cuboid.nX/6, 3*cuboid.nY/4);
+  if (length(q) < 1.1*cuboid.nY/18) {
+    return u_rotate * normalize(make_float2(-q.y, q.x));
+  } else {
+    return make_float2(0);
+  }
+});
+
+auto bulk_mask = materials.mask_of_material(1);
+auto wall_mask = materials.mask_of_material(2);
+auto inflow_mask  = materials.mask_of_material(3);
+auto outflow_mask = materials.mask_of_material(4);
+auto edge_mask = materials.mask_of_material(5);
+
+lattice.apply(Operator(InitializeO(), bulk_mask),
+              Operator(InitializeO(), wall_mask),
+              Operator(InitializeO(), inflow_mask),
+              Operator(InitializeO(), outflow_mask),
+              Operator(InitializeO(), edge_mask));
+cudaDeviceSynchronize();
+
+RenderWindow window("Magnus");
+cudaSurfaceObject_t colormap;
+ColorPalette palette(colormap);
+auto slice = [cuboid] __device__ (int iX, int iY) -> std::size_t {
+               return descriptor::gid(cuboid,iX,cuboid.nY-1-iY);
+             };
+DeviceBuffer<T> moments_rho(cuboid.volume);
+DeviceBuffer<T> moments_u(2*cuboid.volume);
+T* u = moments_u.device();
+std::size_t iStep = 0;
+
+while (window.isOpen()) {
+  lattice.apply(Operator(BgkCollideO(), bulk_mask, tau),
+                Operator(BounceBackFreeSlipO(), wall_mask, WallNormal<0,1>()),
+                Operator(EquilibriumVelocityWallO(), inflow_mask, std::min(iStep*1e-5, 1.)*u_inflow, WallNormal<1,0>()),
+                Operator(EquilibriumDensityWallO(), outflow_mask, 1., WallNormal<-1,0>()),
+                Operator(BounceBackO(), edge_mask));
+  lattice.apply<BouzidiO>(bouzidi.getCount(), bouzidi.getConfig());
+  lattice.stream();
+  if (iStep % 100 == 0) {
+    cudaDeviceSynchronize();
+    lattice.inspect<CollectMomentsF>(bulk_mask, moments_rho.device(), moments_u.device());
+    renderSliceViewToTexture<<<
+      dim3(cuboid.nX / 32 + 1, cuboid.nY / 32 + 1),
+      dim3(32,32)
+    >>>(cuboid.nX, cuboid.nY,
+        slice,
+        [u,u_rotate] __device__ (std::size_t gid) -> float {
+          return length(make_float2(u[2*gid+0], u[2*gid+1])) / u_rotate;
+        },
+        [colormap] __device__ (float x) -> float3 {
+          return colorFromTexture(colormap, clamp(x, 0.f, 1.f));
+        },
+        window.getRenderSurface());
+    window.draw([&]() {
+      ImGui::Begin("Render");
+      palette.interact();
+      ImGui::End();
+    }, [](sf::Event&) { });
+  }
+  ++iStep;
+}
+}
diff --git a/tangle/nozzle.cu b/tangle/nozzle.cu
new file mode 100644
index 0000000..43217dc
--- /dev/null
+++ b/tangle/nozzle.cu
@@ -0,0 +1,73 @@
+#include <LLBM/base.h>
+#include <LLBM/bulk.h>
+#include <LLBM/boundary.h>
+
+#include "util/render_window.h"
+#include "util/texture.h"
+#include "util/colormap.h"
+
+#include "util/volumetric_example.h"
+#include "sampler/velocity_norm.h"
+#include "sampler/curl_norm.h"
+#include "sampler/q_criterion.h"
+
+using T = float;
+using DESCRIPTOR = descriptor::D3Q19;
+
+int main() {
+cudaSetDevice(0);
+
+const descriptor::Cuboid<DESCRIPTOR> cuboid(448, 64, 64);
+Lattice<DESCRIPTOR,T> lattice(cuboid);
+
+CellMaterials<DESCRIPTOR> materials(cuboid, [&cuboid](uint3 p) -> int {
+  if (p.y == 0 || p.y == cuboid.nY-1 || p.z == 0 || p.z == cuboid.nZ-1) {
+    return 2; // boundary cell
+  } else if (p.x == 0) {
+    return 3; // inflow cell
+  } else if (p.x == cuboid.nX-1) {
+    return 4; // outflow cell
+  } else {
+    return 1; // bulk
+  }
+});
+
+auto obstacle = [cuboid] __host__ __device__ (float3 p) -> float {
+                  float3 q = p - make_float3(cuboid.nX/24.2f, cuboid.nY/2, cuboid.nZ/2);
+                  return sdf::ssub(sdf::sphere(make_float2(q.y,q.z), cuboid.nY/T{9}),
+                                   sdf::box(q, make_float3(cuboid.nX/128,cuboid.nY/2,cuboid.nZ/2)),
+                                   5);
+                };
+materials.sdf(obstacle, 0);
+SignedDistanceBoundary bouzidi(lattice, materials, obstacle, 1, 0);
+
+auto bulk_mask     = materials.mask_of_material(1);
+auto boundary_mask = materials.mask_of_material(2);
+auto inflow_mask   = materials.mask_of_material(3);
+auto outflow_mask  = materials.mask_of_material(4);
+
+lattice.apply(Operator(InitializeO(), bulk_mask),
+              Operator(InitializeO(), boundary_mask),
+              Operator(InitializeO(), inflow_mask),
+              Operator(InitializeO(), outflow_mask));
+
+cudaDeviceSynchronize();
+
+VolumetricExample renderer(cuboid);
+renderer.add<CurlNormS>(lattice, bulk_mask, obstacle);
+renderer.add<QCriterionS>(lattice, bulk_mask, obstacle);
+renderer.add<VelocityNormS>(lattice, bulk_mask, obstacle);
+renderer.run([&](std::size_t iStep) {
+  const float tau = 0.501;
+  const float smagorinsky = 0.1;
+  const float inflow = 0.0075;
+  
+  lattice.apply(Operator(SmagorinskyBgkCollideO(), bulk_mask, tau, smagorinsky),
+                Operator(BounceBackO(), boundary_mask),
+                Operator(EquilibriumVelocityWallO(), inflow_mask, std::min(iStep*1e-4, 1.0)*inflow, WallNormal<1,0,0>()),
+                Operator(EquilibriumDensityWallO(), outflow_mask, 1, WallNormal<-1,0,0>()));
+  lattice.apply<BouzidiO>(bouzidi.getCount(), bouzidi.getConfig());
+  
+  lattice.stream();
+});
+}
diff --git a/tangle/sampler/curl_norm.h b/tangle/sampler/curl_norm.h
new file mode 100644
index 0000000..d079a4f
--- /dev/null
+++ b/tangle/sampler/curl_norm.h
@@ -0,0 +1,77 @@
+#pragma once
+
+#include "sampler.h"
+
+#include <LLBM/kernel/collect_moments.h>
+#include <LLBM/kernel/collect_curl.h>
+
+#include <thrust/pair.h>
+#include <thrust/device_vector.h>
+#include <thrust/extrema.h>
+
+template <typename DESCRIPTOR, typename T, typename S, typename SDF>
+class CurlNormS : public Sampler {
+private:
+Lattice<DESCRIPTOR,T,S>& _lattice;
+DeviceBuffer<bool>& _mask;
+SDF _geometry;
+
+DeviceBuffer<float> _moments_rho;
+DeviceBuffer<float> _moments_u;
+DeviceBuffer<float> _curl_norm;
+
+float _scale = 1;
+float _lower = 0;
+float _upper = 1;
+
+public:
+CurlNormS(Lattice<DESCRIPTOR,T,S>& lattice, DeviceBuffer<bool>& mask, SDF geometry):
+  Sampler("Curl norm", lattice.cuboid()),
+  _lattice(lattice),
+  _mask(mask),
+  _geometry(geometry),
+  _moments_rho(lattice.cuboid().volume),
+  _moments_u(DESCRIPTOR::d * lattice.cuboid().volume),
+  _curl_norm(lattice.cuboid().volume)
+{ }
+
+void sample() {
+  _lattice.template inspect<CollectMomentsF>(_mask, _moments_rho.device(), _moments_u.device());
+  _lattice.template inspect<CollectCurlF>(_mask, _moments_u.device(), _sample_surface, _curl_norm.device());
+}
+
+void render(VolumetricRenderConfig& config) {
+  raymarch<<<
+    dim3(config.canvas_size.x / 32 + 1, config.canvas_size.y / 32 + 1),
+    dim3(32, 32)
+  >>>(config,
+      _geometry,
+      [samples=_sample_texture, scale=_scale, lower=_lower, upper=_upper]
+      __device__ (float3 p) -> float {
+        float sample = scale * tex3D<float>(samples, p.x, p.y, p.z);
+        return sample >= lower && sample <= upper ? sample : 0;
+      },
+      [] __device__ (float x) -> float {
+        return x;
+      });
+}
+
+void scale() {
+  auto max = thrust::max_element(thrust::device_pointer_cast(_curl_norm.device()),
+                                 thrust::device_pointer_cast(_curl_norm.device() + _lattice.cuboid().volume));
+  _scale = 1 / max[0];
+}
+
+void interact() {
+  ImGui::SliderFloat("Scale", &_scale, 0.01f, 100.f);
+  ImGui::SameLine();
+  if (ImGui::Button("Auto")) {
+    scale();
+  }
+  ImGui::DragFloatRange2("Bounds", &_lower, &_upper, 0.01f, 0.f, 1.f, "%.2f", "%.2f");
+}
+
+};
+
+template <typename DESCRIPTOR, typename T, typename S, typename SDF>
+CurlNormS(Lattice<DESCRIPTOR,T,S>&, DeviceBuffer<bool>&, SDF) -> CurlNormS<DESCRIPTOR,T,S,SDF>;
diff --git a/tangle/sampler/q_criterion.h b/tangle/sampler/q_criterion.h
new file mode 100644
index 0000000..31eeb61
--- /dev/null
+++ b/tangle/sampler/q_criterion.h
@@ -0,0 +1,90 @@
+#pragma once
+
+#include "sampler.h"
+
+#include <LLBM/kernel/collect_moments.h>
+#include <LLBM/kernel/collect_curl.h>
+#include <LLBM/kernel/collect_q_criterion.h>
+
+#include <thrust/pair.h>
+#include <thrust/device_vector.h>
+#include <thrust/extrema.h>
+
+#include <iostream>
+
+template <typename DESCRIPTOR, typename T, typename S, typename SDF>
+class QCriterionS : public Sampler {
+private:
+Lattice<DESCRIPTOR,T,S>& _lattice;
+DeviceBuffer<bool>& _mask;
+SDF _geometry;
+
+DeviceTexture<float> _curl_buffer;
+cudaTextureObject_t _curl_texture;
+cudaSurfaceObject_t _curl_surface;
+
+DeviceBuffer<float> _moments_rho;
+DeviceBuffer<float> _moments_u;
+DeviceBuffer<float> _curl_norm;
+DeviceBuffer<float> _q;
+
+float _scale = 1;
+float _lower = 0.01;
+float _upper = 1;
+
+public:
+QCriterionS(Lattice<DESCRIPTOR,T,S>& lattice, DeviceBuffer<bool>& mask, SDF geometry):
+  Sampler("Q criterion", lattice.cuboid()),
+  _lattice(lattice),
+  _mask(mask),
+  _geometry(geometry),
+  _curl_buffer(lattice.cuboid()),
+  _curl_texture(_curl_buffer.getTexture()),
+  _curl_surface(_curl_buffer.getSurface()),
+  _moments_rho(lattice.cuboid().volume),
+  _moments_u(DESCRIPTOR::d * lattice.cuboid().volume),
+  _curl_norm(lattice.cuboid().volume),
+  _q(lattice.cuboid().volume)
+{ }
+
+void sample() {
+  _lattice.template inspect<CollectMomentsF>(_mask, _moments_rho.device(), _moments_u.device());
+  _lattice.template inspect<CollectCurlF>(_mask, _moments_u.device(), _curl_surface, _curl_norm.device());
+  _lattice.template inspect<CollectQCriterionF>(_mask, _moments_rho.device(), _moments_u.device(), _curl_norm.device(), _sample_surface, _q.device());
+}
+
+void render(VolumetricRenderConfig& config) {
+  raymarch<<<
+    dim3(config.canvas_size.x / 32 + 1, config.canvas_size.y / 32 + 1),
+    dim3(32, 32)
+  >>>(config,
+      _geometry,
+      [samples=_sample_texture, scale=_scale, lower=_lower, upper=_upper]
+      __device__ (float3 p) -> float {
+        float sample = scale * tex3D<float>(samples, p.x, p.y, p.z);
+        return (sample >= lower) * (sample <= upper) * sample;
+      },
+      [] __device__ (float x) -> float {
+        return (x > 0) * 1;
+      });
+}
+
+void scale() {
+  auto max = thrust::max_element(thrust::device_pointer_cast(_q.device()),
+                                 thrust::device_pointer_cast(_q.device() + _lattice.cuboid().volume));
+  _scale = 1 / max[0];
+}
+
+void interact() {
+  ImGui::SliderFloat("Scale", &_scale, 0.01f, 10000.f);
+  ImGui::SameLine();
+  if (ImGui::Button("Auto")) {
+    scale();
+  }
+  ImGui::DragFloatRange2("Bounds", &_lower, &_upper, 0.01f, 0.01f, 1.f, "%.2f", "%.2f");
+}
+
+};
+
+template <typename DESCRIPTOR, typename T, typename S, typename SDF>
+QCriterionS(Lattice<DESCRIPTOR,T,S>&, DeviceBuffer<bool>&, SDF) -> QCriterionS<DESCRIPTOR,T,S,SDF>;
diff --git a/tangle/sampler/sampler.h b/tangle/sampler/sampler.h
new file mode 100644
index 0000000..ccc50b1
--- /dev/null
+++ b/tangle/sampler/sampler.h
@@ -0,0 +1,32 @@
+#pragma once
+
+#include <LLBM/base.h>
+
+class RenderWindow;
+class VolumetricRenderConfig;
+
+class Sampler {
+protected:
+const std::string _name;
+
+DeviceTexture<float> _sample_buffer;
+cudaTextureObject_t _sample_texture;
+cudaSurfaceObject_t _sample_surface;
+
+public:
+Sampler(std::string name, descriptor::CuboidD<3> cuboid):
+  _name(name),
+  _sample_buffer(cuboid),
+  _sample_texture(_sample_buffer.getTexture()),
+  _sample_surface(_sample_buffer.getSurface())
+ { }
+
+const std::string& getName() const {
+  return _name;
+}
+
+virtual void sample() = 0;
+virtual void render(VolumetricRenderConfig& config) = 0;
+virtual void interact() = 0;
+
+};
diff --git a/tangle/sampler/shear_layer.h b/tangle/sampler/shear_layer.h
new file mode 100644
index 0000000..0d42ee8
--- /dev/null
+++ b/tangle/sampler/shear_layer.h
@@ -0,0 +1,67 @@
+#pragma once
+
+#include "sampler.h"
+
+#include <LLBM/kernel/collect_moments.h>
+#include <LLBM/kernel/collect_shear_layer_normal.h>
+
+template <typename DESCRIPTOR, typename T, typename S, typename SDF>
+class ShearLayerVisibilityS : public Sampler {
+private:
+Lattice<DESCRIPTOR,T,S>& _lattice;
+DeviceBuffer<bool>& _mask;
+SDF _geometry;
+
+DeviceBuffer<float> _moments_rho;
+DeviceBuffer<float> _moments_u;
+DeviceBuffer<float> _shear_normals;
+
+float3 _shear_layer;
+float _lower = 0;
+float _upper = 1;
+bool _center = true;
+
+public:
+ShearLayerVisibilityS(Lattice<DESCRIPTOR,T,S>& lattice, DeviceBuffer<bool>& mask, SDF geometry, float3 shear_layer):
+  Sampler("Shear layer visibility", lattice.cuboid()),
+  _lattice(lattice),
+  _mask(mask),
+  _geometry(geometry),
+  _moments_rho(lattice.cuboid().volume),
+  _moments_u(DESCRIPTOR::d * lattice.cuboid().volume),
+  _shear_normals(DESCRIPTOR::d * lattice.cuboid().volume),
+  _shear_layer(shear_layer)
+{ }
+
+void sample() {
+  _lattice.template inspect<CollectShearLayerNormalsF>(_mask, _moments_rho.device(), _moments_u.device(), _shear_normals.device());
+  _lattice.template helper<CollectShearLayerVisibilityF>(_mask, _shear_normals.device(), _shear_layer, _sample_surface);
+}
+
+void render(VolumetricRenderConfig& config) {
+  raymarch<<<
+    dim3(config.canvas_size.x / 32 + 1, config.canvas_size.y / 32 + 1),
+    dim3(32, 32)
+  >>>(config,
+      _geometry,
+      [samples=_sample_texture, lower=_lower, upper=_upper, center=_center]
+      __device__ (float3 p) -> float {
+        float sample = tex3D<float>(samples, p.x, p.y, p.z);
+        float centered = center ? 0.5 + 0.5*sample : sample;
+        return fabs(sample) >= lower && fabs(sample) <= upper ? fabs(centered) : 0;
+      },
+      [] __device__ (float x) -> float {
+        return x;
+      });
+}
+
+void interact() {
+  ImGui::InputFloat3("Normal", reinterpret_cast<float*>(&_shear_layer));
+  ImGui::Checkbox("Center", &_center);
+  ImGui::DragFloatRange2("Bounds", &_lower, &_upper, 0.01f, 0.f, 1.f, "%.2f", "%.2f");
+}
+
+};
+
+template <typename DESCRIPTOR, typename T, typename S, typename SDF>
+ShearLayerVisibilityS(Lattice<DESCRIPTOR,T,S>&, DeviceBuffer<bool>&, SDF) -> ShearLayerVisibilityS<DESCRIPTOR,T,S,SDF>;
diff --git a/tangle/sampler/velocity_norm.h b/tangle/sampler/velocity_norm.h
new file mode 100644
index 0000000..59ca0a7
--- /dev/null
+++ b/tangle/sampler/velocity_norm.h
@@ -0,0 +1,79 @@
+#pragma once
+
+#include "sampler.h"
+
+#include <LLBM/kernel/collect_moments.h>
+#include <LLBM/kernel/collect_velocity_norm.h>
+
+#include <thrust/pair.h>
+#include <thrust/device_vector.h>
+#include <thrust/extrema.h>
+
+#include <iostream>
+
+template <typename DESCRIPTOR, typename T, typename S, typename SDF>
+class VelocityNormS : public Sampler {
+private:
+Lattice<DESCRIPTOR,T,S>& _lattice;
+DeviceBuffer<bool>& _mask;
+SDF _geometry;
+
+DeviceBuffer<float> _moments_rho;
+DeviceBuffer<float> _moments_u;
+DeviceBuffer<float> _u_norm;
+
+float _scale = 1;
+float _lower = 0;
+float _upper = 1;
+
+public:
+VelocityNormS(Lattice<DESCRIPTOR,T,S>& lattice, DeviceBuffer<bool>& mask, SDF geometry):
+  Sampler("Velocity norm", lattice.cuboid()),
+  _lattice(lattice),
+  _mask(mask),
+  _geometry(geometry),
+  _moments_rho(lattice.cuboid().volume),
+  _moments_u(DESCRIPTOR::d * lattice.cuboid().volume),
+  _u_norm(lattice.cuboid().volume)
+{ }
+
+void sample() {
+  _lattice.template inspect<CollectMomentsF>(_mask, _moments_rho.device(), _moments_u.device());
+  _lattice.template helper<CollectVelocityNormF>(_mask, _moments_u.device(), _sample_surface, _u_norm.device());
+}
+
+void render(VolumetricRenderConfig& config) {
+  raymarch<<<
+    dim3(config.canvas_size.x / 32 + 1, config.canvas_size.y / 32 + 1),
+    dim3(32, 32)
+  >>>(config,
+      _geometry,
+      [samples=_sample_texture, scale=_scale, lower=_lower, upper=_upper]
+      __device__ (float3 p) -> float {
+        float sample = scale * tex3D<float>(samples, p.x, p.y, p.z);
+        return sample >= lower && sample <= upper ? sample : 0;
+      },
+      [] __device__ (float x) -> float {
+        return x;
+      });
+}
+
+void scale() {
+  auto max = thrust::max_element(thrust::device_pointer_cast(_u_norm.device()),
+                                 thrust::device_pointer_cast(_u_norm.device() + _lattice.cuboid().volume));
+  _scale = 1 / max[0];
+}
+
+void interact() {
+  ImGui::SliderFloat("Scale", &_scale, 0.01f, 100.f);
+  ImGui::SameLine();
+  if (ImGui::Button("Auto")) {
+    scale();
+  }
+  ImGui::DragFloatRange2("Bounds", &_lower, &_upper, 0.01f, 0.f, 1.f, "%.2f", "%.2f");
+}
+
+};
+
+template <typename DESCRIPTOR, typename T, typename S, typename SDF>
+VelocityNormS(Lattice<DESCRIPTOR,T,S>&, DeviceBuffer<bool>&, SDF) -> VelocityNormS<DESCRIPTOR,T,S,SDF>;
diff --git a/tangle/taylor-couette.cu b/tangle/taylor-couette.cu
new file mode 100644
index 0000000..86f739e
--- /dev/null
+++ b/tangle/taylor-couette.cu
@@ -0,0 +1,75 @@
+#include <LLBM/base.h>
+#include <LLBM/bulk.h>
+#include <LLBM/boundary.h>
+
+#include "util/render_window.h"
+#include "util/texture.h"
+#include "util/colormap.h"
+
+#include "util/volumetric_example.h"
+#include "sampler/velocity_norm.h"
+#include "sampler/curl_norm.h"
+#include "sampler/shear_layer.h"
+
+using T = float;
+using DESCRIPTOR = descriptor::D3Q19;
+
+int main() {
+cudaSetDevice(0);
+
+const descriptor::Cuboid<DESCRIPTOR> cuboid(320, 96, 96);
+Lattice<DESCRIPTOR,T> lattice(cuboid);
+
+CellMaterials<DESCRIPTOR> materials(cuboid, [&cuboid](uint3 p) -> int {
+  if (p.x == 0 || p.x == cuboid.nX-1) {
+    return 2;
+  } else {
+    return 1;
+  }
+});
+
+auto inner_cylinder = [cuboid] __host__ __device__ (float3 p) -> float {
+                        float3 q = p - make_float3(0, cuboid.nY/2, cuboid.nZ/2);
+                        return sdf::sphere(make_float2(q.y,q.z), cuboid.nY/T{4.5});
+                      };
+auto geometry = [cuboid,inner_cylinder] __host__ __device__ (float3 p) -> float {
+                  float3 q = p - make_float3(0, cuboid.nY/2, cuboid.nZ/2);
+                  return sdf::add(-sdf::sphere(make_float2(q.y,q.z), cuboid.nY/T{2.14}), inner_cylinder(p));
+                };
+materials.sdf(geometry, 0);
+SignedDistanceBoundary bouzidi(lattice, materials, geometry, 1, 0);
+
+const float wall = 0.2;
+
+bouzidi.setVelocity([cuboid,wall](float3 p) -> float3 {
+  float3 q = p - make_float3(0, cuboid.nY/2, cuboid.nZ/2);
+  if (length(make_float2(q.y,q.z)) < cuboid.nY/T{2.5}) {
+    return wall * normalize(make_float3(0, -q.z, q.y));
+  } else {
+    return make_float3(0);
+  }
+});
+
+auto bulk_mask = materials.mask_of_material(1);
+auto bulk_list = materials.list_of_material(1);
+auto wall_mask = materials.mask_of_material(2);
+auto wall_list = materials.list_of_material(2);
+
+lattice.apply<InitializeO>(bulk_list);
+lattice.apply<InitializeO>(wall_list);
+
+cudaDeviceSynchronize();
+
+VolumetricExample renderer(cuboid);
+renderer.add<VelocityNormS>(lattice, bulk_mask, inner_cylinder);
+renderer.add<ShearLayerVisibilityS>(lattice, bulk_mask, inner_cylinder, make_float3(1,0,0));
+renderer.run([&](std::size_t iStep) {
+  const float tau = 0.55;
+  
+  lattice.apply<BgkCollideO>(bulk_list, tau);
+  lattice.apply<BounceBackO>(wall_list);
+  lattice.apply<BouzidiO>(bouzidi.getCount(), bouzidi.getConfig());
+  
+  lattice.stream();
+});
+}
diff --git a/tangle/tmp/noise_overview.png b/tangle/tmp/noise_overview.png
new file mode 100644
index 0000000..2e657ac
--- /dev/null
+++ b/tangle/tmp/noise_overview.png
diff --git a/tangle/tmp/test_noise.png b/tangle/tmp/test_noise.png
new file mode 100644
index 0000000..e9056f7
--- /dev/null
+++ b/tangle/tmp/test_noise.png
diff --git a/tangle/util/camera.h b/tangle/util/camera.h
new file mode 100644
index 0000000..0a793d3
--- /dev/null
+++ b/tangle/util/camera.h
@@ -0,0 +1,86 @@
+#include <cuda-samples/Common/helper_math.h>
+#include "SFML/Window/Event.hpp"
+
+class Camera {
+private:
+  float _distance;
+  float _phi;
+  float _psi;
+  float3 _target;
+  float3 _eye;
+  float3 _direction;
+  bool _dragging;
+  bool _moving;
+  float2 _lastMouse;
+
+public:
+  Camera(float3 target, float distance, float phi, float psi):
+    _distance(distance),
+    _phi(phi),
+    _psi(psi),
+    _target(target),
+    _dragging(false),
+    _moving(false) {
+    update();
+  }
+
+  void update() {
+    _eye = _target + make_float3(_distance*sin(_psi)*cos(_phi), _distance*sin(_psi)*sin(_phi), _distance*cos(_psi));
+    _direction = normalize(_target - _eye);
+  }
+
+  void handle(sf::Event& event) {
+    switch (event.type) {
+    case sf::Event::MouseWheelMoved:
+      _distance -= event.mouseWheel.delta * 10;
+      break;
+    case sf::Event::MouseButtonPressed:
+      if (event.mouseButton.button == sf::Mouse::Left) {
+        _dragging = true;
+        _lastMouse = make_float2(event.mouseButton.x, event.mouseButton.y);
+      } else if (event.mouseButton.button == sf::Mouse::Right) {
+        _moving = true;
+        _lastMouse = make_float2(event.mouseButton.x, event.mouseButton.y);
+      }
+      break;
+    case sf::Event::MouseButtonReleased:
+      if (event.mouseButton.button == sf::Mouse::Left) {
+        _dragging = false;
+      } else if (event.mouseButton.button == sf::Mouse::Right) {
+        _moving = false;
+      }
+      break;
+    case sf::Event::MouseMoved:
+      if (_dragging) {
+        float2 mouse = make_float2(event.mouseMove.x, event.mouseMove.y);
+        float2 delta = mouse - _lastMouse;
+        _lastMouse = mouse;
+        _phi += 0.4*delta.x * 2*M_PI/360;
+        if (delta.y > 0 && _psi <= M_PI-2*M_PI/60) {
+          _psi += 0.4*delta.y * M_PI/180;
+        } else if (delta.y < 0 && _psi >= 2*M_PI/60) {
+          _psi += 0.4*delta.y * M_PI/180;
+        }
+      }
+      if (_moving) {
+        float2 mouse = make_float2(event.mouseMove.x, event.mouseMove.y);
+        float2 delta = mouse - _lastMouse;
+        _lastMouse = mouse;
+        float3 forward = normalize(_target - _eye);
+    	   float3 right   = normalize(cross(make_float3(0.f, 0.f, -1.f), forward));
+        float3 up      = cross(right, forward);
+        _target += 0.4*right*delta.x - 0.4*up*delta.y;
+      }
+      break;
+    }
+    update();
+  }
+
+  float3 getDirection() const {
+    return _direction;
+  }
+
+  float3 getEyePosition() const {
+    return _eye;
+  }
+};
diff --git a/tangle/util/colormap.h b/tangle/util/colormap.h
new file mode 100644
index 0000000..c01cea8
--- /dev/null
+++ b/tangle/util/colormap.h
@@ -0,0 +1,38 @@
+#pragma once
+#include "assets.h"
+#include "texture.h"
+
+#include <imgui.h>
+#include <imgui-SFML.h>
+#include <SFML/Graphics.hpp>
+
+struct ColorPalette {
+  const assets::File* current;
+  sf::Texture texture;
+
+  ColorPalette(cudaSurfaceObject_t& palette) {
+    current = &assets::palette::files[5];
+    texture.loadFromMemory(current->data, current->size);
+    palette = bindTextureToCuda(texture);
+  }
+
+  void interact();
+};
+
+void ColorPalette::interact() {
+  if (ImGui::BeginCombo("Color palette", current->name.c_str())) {
+    for (unsigned i=0; i < assets::palette::file_count; ++i) {
+      bool is_selected = (current == &assets::palette::files[i]);
+      if (ImGui::Selectable(assets::palette::files[i].name.c_str(), is_selected)) {
+        current = &assets::palette::files[i];
+        texture.loadFromMemory(current->data, current->size);
+        break;
+      }
+      if (is_selected) {
+        ImGui::SetItemDefaultFocus();
+      }
+    }
+    ImGui::EndCombo();
+  }
+  ImGui::Image(texture, sf::Vector2f(400.,40.));
+}
diff --git a/tangle/util/noise.h b/tangle/util/noise.h
new file mode 100644
index 0000000..22ad548
--- /dev/null
+++ b/tangle/util/noise.h
@@ -0,0 +1,39 @@
+#pragma once
+#include "assets.h"
+#include "texture.h"
+
+#include <imgui.h>
+#include <imgui-SFML.h>
+#include <SFML/Graphics.hpp>
+
+struct NoiseSource {
+  const assets::File* current;
+  sf::Texture texture;
+
+  NoiseSource(cudaSurfaceObject_t& noise) {
+    current = &assets::noise::files[0];
+    texture.loadFromMemory(current->data, current->size);
+    noise = bindTextureToCuda(texture);
+  }
+
+  void interact();
+};
+
+void NoiseSource::interact() {
+  ImGui::Image(texture, sf::Vector2f(32,20));
+  ImGui::SameLine();
+  if (ImGui::BeginCombo("Noise", current->name.c_str())) {
+    for (unsigned i=0; i < assets::noise::file_count; ++i) {
+      bool is_selected = (current == &assets::noise::files[i]);
+      if (ImGui::Selectable(assets::noise::files[i].name.c_str(), is_selected)) {
+        current = &assets::noise::files[i];
+        texture.loadFromMemory(current->data, current->size);
+        break;
+      }
+      if (is_selected) {
+        ImGui::SetItemDefaultFocus();
+      }
+    }
+    ImGui::EndCombo();
+  }
+}
diff --git a/tangle/util/render_window.h b/tangle/util/render_window.h
new file mode 100644
index 0000000..4392150
--- /dev/null
+++ b/tangle/util/render_window.h
@@ -0,0 +1,95 @@
+#pragma once
+
+#include <SFML/Graphics.hpp>
+#include <SFML/Graphics/Image.hpp>
+
+#include <imgui.h>
+#include <imgui-SFML.h>
+
+#include "texture.h"
+#include "assets.h"
+
+class RenderWindow {
+private:
+sf::RenderWindow _window;
+
+sf::Sprite          _render_sprite;
+sf::Texture         _render_texture;
+cudaSurfaceObject_t _render_surface;
+sf::Rect<int>       _render_texture_view;
+
+sf::Shader _blur_shader;
+bool _blur = false;
+
+sf::Clock _ui_delta_clock;
+
+public:
+RenderWindow(std::string name):
+  _window(sf::VideoMode(800, 600), name) {
+  _render_texture.create(sf::VideoMode::getDesktopMode().width, sf::VideoMode::getDesktopMode().height);
+  _render_surface = bindTextureToCuda(_render_texture);
+  _render_sprite.setTexture(_render_texture);
+  _render_texture_view = sf::Rect<int>(0,0,_window.getSize().x,_window.getSize().y);
+  _render_sprite.setTextureRect(_render_texture_view);
+  _window.setView(sf::View(sf::Vector2f(_render_texture_view.width/2, _render_texture_view.height/2),
+                           sf::Vector2f(_window.getSize().x, _window.getSize().y)));
+  _window.setVerticalSyncEnabled(true);
+  _blur_shader.loadFromMemory(std::string(reinterpret_cast<const char*>(assets::shader::file_blur_frag)), sf::Shader::Fragment);
+  _blur_shader.setUniform("texture", sf::Shader::CurrentTexture);
+  ImGui::SFML::Init(_window);
+  ImGuiIO& io = ImGui::GetIO();
+  io.MouseDrawCursor = true;
+};
+
+bool isOpen() const {
+  return _window.isOpen();
+}
+
+void setBlur(bool state) {
+  _blur = state;
+}
+
+template <typename UI, typename MOUSE>
+void draw(UI ui, MOUSE mouse);
+
+cudaSurfaceObject_t getRenderSurface() {
+  return _render_surface;
+}
+
+sf::Rect<int> getRenderView() {
+  return _render_texture_view;
+}
+
+};
+
+template <typename UI, typename MOUSE>
+void RenderWindow::draw(UI ui, MOUSE mouse) {
+  sf::Event event;
+  while (_window.pollEvent(event)) {
+    ImGui::SFML::ProcessEvent(event);
+    if (event.type == sf::Event::Closed) {
+      _window.close();
+    }
+    if (event.type == sf::Event::Resized) {
+      _render_texture_view = sf::Rect<int>(0,0,event.size.width,event.size.height);
+      _render_sprite.setTextureRect(_render_texture_view);
+      sf::View view(sf::Vector2f(_render_texture_view.width/2, _render_texture_view.height/2),
+                    sf::Vector2f(event.size.width, event.size.height));
+      _window.setView(view);
+    }
+    if (!ImGui::GetIO().WantCaptureMouse) {
+      mouse(event);
+    }
+  }
+
+  ImGui::SFML::Update(_window, _ui_delta_clock.restart());
+  ui();
+  _window.clear();
+  if (_blur) {
+    _window.draw(_render_sprite, &_blur_shader);
+  } else {
+    _window.draw(_render_sprite);
+  }
+  ImGui::SFML::Render(_window);
+  _window.display();
+}
diff --git a/tangle/util/texture.h b/tangle/util/texture.h
new file mode 100644
index 0000000..f2e0440
--- /dev/null
+++ b/tangle/util/texture.h
@@ -0,0 +1,25 @@
+#pragma once
+
+#include <cstring>
+#include <SFML/Graphics.hpp>
+#include <cuda_gl_interop.h>
+#include <LLBM/memory.h>
+
+cudaSurfaceObject_t bindTextureToCuda(sf::Texture& texture) {
+  GLuint gl_tex_handle = texture.getNativeHandle();
+  cudaGraphicsResource* cuda_tex_handle;
+  cudaArray* buffer;
+
+  cudaGraphicsGLRegisterImage(&cuda_tex_handle, gl_tex_handle, GL_TEXTURE_2D, cudaGraphicsRegisterFlagsNone);
+  cudaGraphicsMapResources(1, &cuda_tex_handle, 0);
+  cudaGraphicsSubResourceGetMappedArray(&buffer, cuda_tex_handle, 0, 0);
+
+  cudaResourceDesc resDesc;
+  resDesc.resType = cudaResourceTypeArray;
+
+  resDesc.res.array.array = buffer;
+  cudaSurfaceObject_t cudaSurfaceObject = 0;
+  cudaCreateSurfaceObject(&cudaSurfaceObject, &resDesc);
+
+  return cudaSurfaceObject;
+}
diff --git a/tangle/util/timer.h b/tangle/util/timer.h
new file mode 100644
index 0000000..fd5d832
--- /dev/null
+++ b/tangle/util/timer.h
@@ -0,0 +1,19 @@
+#pragma once
+#include <chrono>
+
+namespace timer {
+
+std::chrono::time_point<std::chrono::steady_clock> now() {
+  return std::chrono::steady_clock::now();
+}
+
+double secondsSince(
+  std::chrono::time_point<std::chrono::steady_clock>& pit) {
+  return std::chrono::duration_cast<std::chrono::duration<double>>(now() - pit).count();
+}
+
+double mlups(std::size_t nCells, std::size_t nSteps, std::chrono::time_point<std::chrono::steady_clock>& start) {
+  return nCells * nSteps / (secondsSince(start) * 1e6);
+}
+
+}
diff --git a/tangle/util/volumetric_example.h b/tangle/util/volumetric_example.h
new file mode 100644
index 0000000..cc6a24e
--- /dev/null
+++ b/tangle/util/volumetric_example.h
@@ -0,0 +1,121 @@
+#pragma once
+#include <LLBM/volumetric.h>
+
+#include "camera.h"
+#include "texture.h"
+#include "colormap.h"
+#include "noise.h"
+#include "render_window.h"
+#include "../sampler/sampler.h"
+
+class VolumetricExample : public RenderWindow {
+private:
+std::vector<std::unique_ptr<Sampler>> _sampler;
+Sampler* _current = nullptr;
+
+Camera _camera;
+VolumetricRenderConfig _config;
+ColorPalette _palette;
+NoiseSource _noise;
+
+int _steps_per_second = 100;
+int _samples_per_second = 30;
+
+public:
+VolumetricExample(descriptor::CuboidD<3> cuboid):
+  RenderWindow("LiterateLB"),
+  _camera(make_float3(cuboid.nX/2,cuboid.nY/2,cuboid.nZ/2), cuboid.nX, M_PI/2, M_PI/2),
+  _config(cuboid),
+  _palette(_config.palette),
+  _noise(_config.noise)
+{
+  _config.canvas = this->getRenderSurface();
+  this->setBlur(_config.apply_blur);
+}
+
+template <template<typename...> class SAMPLER, typename... ARGS>
+void add(ARGS&&... args) {
+  _sampler.emplace_back(new SAMPLER(std::forward<ARGS>(args)...));
+  _current = _sampler.back().get();
+}
+
+template <typename TIMESTEP>
+void run(TIMESTEP step) {
+  sf::Clock last_sample;
+  sf::Clock last_frame;
+  std::size_t iStep = 0;
+  volatile bool simulate = true;
+
+  sf::Thread simulation([&]() {
+    while (this->isOpen()) {
+      if (last_sample.getElapsedTime().asSeconds() > 1.0 / _samples_per_second) {
+        _current->sample();
+        cudaStreamSynchronize(cudaStreamPerThread);
+        last_sample.restart();
+        if (simulate) {
+          for (unsigned i=0; i < (1.0 / _samples_per_second) * _steps_per_second; ++i) {
+            step(iStep++);
+          }
+        }
+      }
+    }
+  });
+  simulation.launch();
+
+  while (this->isOpen()) {
+    this->draw(
+      [&](){
+        ImGui::Begin("Simulation", 0, ImGuiWindowFlags_AlwaysAutoResize);
+        if (ImGui::BeginCombo("Source", _current->getName().c_str())) {
+          for (auto& option : _sampler) {
+            if (ImGui::Selectable(option->getName().c_str(), _current == option.get())) {
+              _current = option.get();
+            }
+          }
+          ImGui::EndCombo();
+        }
+        _current->interact();
+        ImGui::SliderInt("Timestep/s", &_steps_per_second, 1, 2000);
+        ImGui::SliderInt("Samples/s", &_samples_per_second, 1, 60);
+        if (simulate) {
+          simulate = !ImGui::Button("Pause");
+        } else {
+          simulate =  ImGui::Button("Continue");
+        }
+        ImGui::End();
+        ImGui::Begin("Render", 0, ImGuiWindowFlags_AlwaysAutoResize);
+        ImGui::SliderFloat("Brightness", &_config.brightness, 0.1f, 2.f);
+        ImGui::SliderFloat("Delta", &_config.delta, 0.05f, 2.f);
+        ImGui::SliderFloat("Transparency", &_config.transparency, 0.001f, 1.f);
+        _palette.interact();
+        if (ImGui::CollapsingHeader("Details")) {
+          ImGui::Checkbox("Align slices to view", &_config.align_slices_to_view);
+          ImGui::SameLine();
+          ImGui::Checkbox("Jitter", &_config.apply_noise);
+          ImGui::SameLine();
+          ImGui::Checkbox("Blur", &_config.apply_blur);
+          this->setBlur(_config.apply_blur);
+          if (_config.apply_noise) {
+            _noise.interact();
+          }
+        }
+        ImGui::End();
+      },
+      [&](sf::Event& event) {
+        _camera.handle(event);
+        _config.eye_pos = _camera.getEyePosition();
+        _config.eye_dir = _camera.getDirection();
+        _config.canvas_size = make_uint2(this->getRenderView().width, this->getRenderView().height);
+      }
+    );
+    if (last_frame.getElapsedTime().asSeconds() > 1.0 / _samples_per_second) {
+      _current->render(_config);
+      cudaStreamSynchronize(cudaStreamPerThread);
+      last_frame.restart();
+    }
+  }
+
+  simulation.wait();
+}
+
+};