Small changes to channel example

3 files changed, 62 insertions, 39 deletions

diff --git a/lbm.org b/lbm.org
index 168427b..b21ca1d 100644
--- a/lbm.org
+++ b/lbm.org
@@ -157,7 +157,7 @@ speaking compiling and testing the examples requires neither Emacs nor Python.
 
 Note that I started developing this as a beginner in both Emacs and Org mode so some aspects of this document
 may be more clunky than necessary. Most of the complexity stems from maintaining a Python session for the
-generation of optimized GPU kernel functions using the SymPy CAS:
+generation of optimized GPU kernel functions using the SymPy CAS.
 * Lattice
 The cartesian grids used for the spatial discretization are commonly described as =DXQY= lattices where =X= is the number of spatial dimensions
 and =Y= is the number of discrete velocities $\xi_i$. e.g. a  =D2Q9= lattice is two dimensional and stores nine population values per cell. Each population has
@@ -173,7 +173,6 @@ class Descriptor:
         self.c = [ Matrix(eval(row[0]))    for row in data ]
         self.w = [ Rational(f.numerator, f.denominator)
                    for f in [ Fraction(row[1]) for row in data ] ]
-        
         self.d = self.c[0].shape[0]
         self.q = len(self.c)
         self.c_s = sqrt(Rational(1,3))
@@ -261,14 +260,14 @@ class CodeBlockPrinter(C11CodePrinter):
         self.custom_assignment = custom_assignment
         for f in custom_functions:
             self._kf[f] = f
-    
+
     def _print_Indexed(self, expr):
         assert len(expr.indices) == 1
         if expr.base.name[0] == 'f':
             return f"{expr.base.name}[{expr.indices[0]}]"
         else:
             return f"{expr.base.name}_{expr.indices[0]}"
-    
+
     def _print_Float(self, flt):
         return "T{%s}" % str(flt.evalf())
 
@@ -1666,6 +1665,10 @@ A convenient way for generating SDFs for arbitrary shapes is to construct them b
 various primitives such as spheres and boxes using boolean operators such as addition and
 substraction.
 
+A comprensive listing of different such functions is available on [[https://iquilezles.org/www/articles/distfunctions/distfunctions.htm][iquilezles.org]].
+The remainder of this section translates some of these shader functions into C++
+functions executable both on the host and the device.
+
 #+BEGIN_SRC cpp :tangle tangle/LLBM/sdf.h
 namespace sdf {
 #+END_SRC
@@ -1748,6 +1751,17 @@ __device__ __host__ float sintersect(float a, float b, float k) {
 #+END_SRC
 
 #+BEGIN_SRC cpp :tangle tangle/LLBM/sdf.h
+__device__ __host__ float3 twisted(float3 p, float k) {
+  float c = cos(k*p.y);
+  float s = sin(k*p.y);
+  float3  q = make_float3(0,0,p.y);
+  q.x = p.x*c + p.z*-s;
+  q.y = p.x*s + p.z* c;
+  return q;
+}
+#+END_SRC
+
+#+BEGIN_SRC cpp :tangle tangle/LLBM/sdf.h
 }
 #+END_SRC
 
@@ -3086,20 +3100,15 @@ def readColormap(src):
     for s in root.findall('.//Point'):
         values.append(float(s.attrib['x']))
         colors.append((float(s.attrib['r']), float(s.attrib['g']), float(s.attrib['b'])))
-    
     colormap = [ ]
     if values[0] != 0:
         colormap.append((0, colors[0]))
-    
     for pos, color in zip(values, colors):
         colormap.append((pos, color))
-    
     if values[-1] != 1:
         colormap.append((1, colors[-1]))
-    
     return colormap
 
-
 def renderColormap(cmap, path, reversed = False):
     gradient = np.linspace(1, 0, 1000) if reversed else np.linspace(0, 1, 1000)
     gradient = np.vstack((gradient, gradient))
@@ -5377,10 +5386,10 @@ magnus effect and the formation of a Kármán vortex street.
 <video style="width:100%" src="https://literatelb.org/media/magnus.webm" playsinline muted controls loop/>
 #+END_EXPORT
 
-** Flow around a Sphere
+** Channel with obstacle
 This example models a channel flow around a spherical obstacle.
 
-#+BEGIN_SRC cpp :tangle tangle/channel-with-sphere.cu
+#+BEGIN_SRC cpp :tangle tangle/channel.cu
 <<example-headers>>
 
 #include "util/volumetric_example.h"
@@ -5401,15 +5410,15 @@ auto current = cuda::device::current::get();
 
 After including the relevant headers we construct the D3Q19 lattice.
 
-#+BEGIN_SRC cpp :tangle tangle/channel-with-sphere.cu
-const descriptor::Cuboid<DESCRIPTOR> cuboid(448, 64, 64);
+#+BEGIN_SRC cpp :tangle tangle/channel.cu
+const descriptor::Cuboid<DESCRIPTOR> cuboid(500, 100, 100);
 Lattice<DESCRIPTOR,T> lattice(cuboid);
 #+END_SRC
 
 Before constructing the sphere we set up a basic square channel with
 free slip walls and separate in- and outlets.
 
-#+BEGIN_SRC cpp :tangle tangle/channel-with-sphere.cu
+#+BEGIN_SRC cpp :tangle tangle/channel.cu
 CellMaterials<DESCRIPTOR> materials(cuboid, [&cuboid](uint3 p) -> int {
   if (p.z == 0 || p.z == cuboid.nZ-1) {
     return 2; // boundary cell
@@ -5428,7 +5437,7 @@ CellMaterials<DESCRIPTOR> materials(cuboid, [&cuboid](uint3 p) -> int {
 As the free slip bounce back condition is not well defined at the edges we use
 plain bounce back instead.
 
-#+BEGIN_SRC cpp :tangle tangle/channel-with-sphere.cu
+#+BEGIN_SRC cpp :tangle tangle/channel.cu
 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);
@@ -5437,22 +5446,23 @@ for (std::size_t iX=0; iX < cuboid.nX; ++iX) {
 }
 #+END_SRC
 
-The sphere is modelled using a signed distance function from which we
+The obstacle is modelled using a signed distance function from which we
 generate a configuraton of the interpolated bounce back condition.
 
-#+BEGIN_SRC cpp :tangle tangle/channel-with-sphere.cu
+#+BEGIN_SRC cpp :tangle tangle/channel.cu
 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});
-                };
+  p -= make_float3(cuboid.nX/5, cuboid.nY/2, cuboid.nZ/2);
+  float3 q = sdf::twisted(p, 0.01);
+  return sdf::sphere(p, cuboid.nY/3.5) + sin(0.2*q.x)*sin(0.2*q.y)*sin(0.2*q.z);
+};
+
 materials.sdf(obstacle, 0);
 SignedDistanceBoundary bouzidi(lattice, materials, obstacle, 1, 0);
 #+END_SRC
 
-At this point themasks for operator application can be generated followed by
-initializing the lattice data.
+At this point the masks for operator application can be generated.
 
-#+BEGIN_SRC cpp :tangle tangle/channel-with-sphere.cu
+#+BEGIN_SRC cpp :tangle tangle/channel.cu
 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);
@@ -5463,16 +5473,17 @@ auto edge_mask    = materials.mask_of_material(6);
 cuda::synchronize(current);
 #+END_SRC
 
-In order to model the flow we employ standard BGK collisions in the bulk. The channel
-walls are modelled using (free slip) bounce back and the SDF-described sphere is
-represented with interpolated bounce back.
+In order to model the flow we employ Smagorinsky BGK collisions in the bulk. The
+channel walls are modelled using (free slip) bounce back and the SDF-described
+obstacle is represented with interpolated bounce back.
 
-#+NAME: channel-with-sphere-simulation-step
+#+NAME: channel-simulation-step
 #+BEGIN_SRC cpp :eval no
-const float tau = 0.51;
-const float inflow = 0.05;
+const float tau = 0.501;
+const float smagorinsky = 0.1;
+const float inflow = 0.04;
 
-lattice.apply(Operator(BgkCollideO(), bulk_mask, tau),
+lattice.apply(Operator(SmagorinskyBgkCollideO(), bulk_mask, tau, smagorinsky),
               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>()),
@@ -5486,13 +5497,13 @@ lattice.stream();
 Finally the volumetric renderer is used to control the simulation
 loop and to provide velocity, curl and Q criterion visualizations.
 
-#+BEGIN_SRC cpp :tangle tangle/channel-with-sphere.cu
+#+BEGIN_SRC cpp :tangle tangle/channel.cu
 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) {
-  <<channel-with-sphere-simulation-step>>
+  <<channel-simulation-step>>
 });
 }
 #+END_SRC
diff --git a/tangle/LLBM/sdf.h b/tangle/LLBM/sdf.h
index 656c109..52db8cf 100644
--- a/tangle/LLBM/sdf.h
+++ b/tangle/LLBM/sdf.h
@@ -63,6 +63,15 @@ __device__ __host__ float sintersect(float a, float b, float k) {
   return lerp(b, a, h) + k*h*(1.f-h);
 }
 
+__device__ __host__ float3 twisted(float3 p, float k) {
+  float c = cos(k*p.y);
+  float s = sin(k*p.y);
+  float3  q = make_float3(0,0,p.y);
+  q.x = p.x*c + p.z*-s;
+  q.y = p.x*s + p.z* c;
+  return q;
+}
+
 }
 
 template <typename SDF>
diff --git a/tangle/channel-with-sphere.cu b/tangle/channel.cu
index 280a142..6a0f072 100644
--- a/tangle/channel-with-sphere.cu
+++ b/tangle/channel.cu
@@ -21,7 +21,7 @@ if (cuda::device::count() == 0) {
 }
 auto current = cuda::device::current::get();
 
-const descriptor::Cuboid<DESCRIPTOR> cuboid(448, 64, 64);
+const descriptor::Cuboid<DESCRIPTOR> cuboid(500, 100, 100);
 Lattice<DESCRIPTOR,T> lattice(cuboid);
 
 CellMaterials<DESCRIPTOR> materials(cuboid, [&cuboid](uint3 p) -> int {
@@ -46,9 +46,11 @@ for (std::size_t iX=0; iX < cuboid.nX; ++iX) {
 }
 
 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});
-                };
+  p -= make_float3(cuboid.nX/5, cuboid.nY/2, cuboid.nZ/2);
+  float3 q = sdf::twisted(p, 0.01);
+  return sdf::sphere(p, cuboid.nY/3.5) + sin(0.2*q.x)*sin(0.2*q.y)*sin(0.2*q.z);
+};
+
 materials.sdf(obstacle, 0);
 SignedDistanceBoundary bouzidi(lattice, materials, obstacle, 1, 0);
 
@@ -66,10 +68,11 @@ 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.05;
+  const float tau = 0.501;
+  const float smagorinsky = 0.1;
+  const float inflow = 0.04;
   
-  lattice.apply(Operator(BgkCollideO(), bulk_mask, tau),
+  lattice.apply(Operator(SmagorinskyBgkCollideO(), bulk_mask, tau, smagorinsky),
                 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>()),