From 186bfae321f9d34b702a7c31b420eed4d10e7a50 Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Sun, 19 Sep 2021 21:02:52 +0200
Subject: Small changes to channel example
---
lbm.org | 75 ++++++++++++++++++++++----------------
tangle/LLBM/sdf.h | 9 +++++
tangle/channel-with-sphere.cu | 82 -----------------------------------------
tangle/channel.cu | 85 +++++++++++++++++++++++++++++++++++++++++++
4 files changed, 137 insertions(+), 114 deletions(-)
delete mode 100644 tangle/channel-with-sphere.cu
create mode 100644 tangle/channel.cu
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
@@ -1747,6 +1750,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.
#+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
<>
#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 cuboid(448, 64, 64);
+#+BEGIN_SRC cpp :tangle tangle/channel.cu
+const descriptor::Cuboid cuboid(500, 100, 100);
Lattice 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 materials(cuboid, [&cuboid](uint3 p) -> int {
if (p.z == 0 || p.z == cuboid.nZ-1) {
return 2; // boundary cell
@@ -5428,7 +5437,7 @@ CellMaterials 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(lattice, bulk_mask, obstacle);
renderer.add(lattice, bulk_mask, obstacle);
renderer.add(lattice, bulk_mask, obstacle);
renderer.run([&](std::size_t iStep) {
- <>
+ <>
});
}
#+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
diff --git a/tangle/channel-with-sphere.cu b/tangle/channel-with-sphere.cu
deleted file mode 100644
index 280a142..0000000
--- a/tangle/channel-with-sphere.cu
+++ /dev/null
@@ -1,82 +0,0 @@
-#include
-#include
-#include
-
-#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() {
-if (cuda::device::count() == 0) {
- std::cerr << "No CUDA devices on this system" << std::endl;
- return -1;
-}
-auto current = cuda::device::current::get();
-
-const descriptor::Cuboid cuboid(448, 64, 64);
-Lattice lattice(cuboid);
-
-CellMaterials 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);
-
-cuda::synchronize(current);
-
-VolumetricExample renderer(cuboid);
-renderer.add(lattice, bulk_mask, obstacle);
-renderer.add(lattice, bulk_mask, obstacle);
-renderer.add(lattice, bulk_mask, obstacle);
-renderer.run([&](std::size_t iStep) {
- const float tau = 0.51;
- const float inflow = 0.05;
-
- 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(bouzidi.getCount(), bouzidi.getConfig());
-
- lattice.stream();
-});
-}
diff --git a/tangle/channel.cu b/tangle/channel.cu
new file mode 100644
index 0000000..6a0f072
--- /dev/null
+++ b/tangle/channel.cu
@@ -0,0 +1,85 @@
+#include
+#include
+#include
+
+#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() {
+if (cuda::device::count() == 0) {
+ std::cerr << "No CUDA devices on this system" << std::endl;
+ return -1;
+}
+auto current = cuda::device::current::get();
+
+const descriptor::Cuboid cuboid(500, 100, 100);
+Lattice lattice(cuboid);
+
+CellMaterials 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 {
+ 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);
+
+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);
+
+cuda::synchronize(current);
+
+VolumetricExample renderer(cuboid);
+renderer.add(lattice, bulk_mask, obstacle);
+renderer.add(lattice, bulk_mask, obstacle);
+renderer.add(lattice, bulk_mask, obstacle);
+renderer.run([&](std::size_t iStep) {
+ const float tau = 0.501;
+ const float smagorinsky = 0.1;
+ const float inflow = 0.04;
+
+ 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>()),
+ Operator(EquilibriumDensityWallO(), outflow_mask, 1, WallNormal<-1,0,0>()),
+ Operator(BounceBackO(), edge_mask));
+ lattice.apply(bouzidi.getCount(), bouzidi.getConfig());
+
+ lattice.stream();
+});
+}
--
cgit v1.2.3