summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAdrian Kummerlaender2021-09-19 21:02:52 +0200
committerAdrian Kummerlaender2021-09-19 21:02:52 +0200
commit186bfae321f9d34b702a7c31b420eed4d10e7a50 (patch)
tree4e8e6a691193dd0ce9555d2e29f4414336d17cda
parent32dd41a728ce10113032e20955ba08f8de449857 (diff)
downloadLiterateLB-186bfae321f9d34b702a7c31b420eed4d10e7a50.tar
LiterateLB-186bfae321f9d34b702a7c31b420eed4d10e7a50.tar.gz
LiterateLB-186bfae321f9d34b702a7c31b420eed4d10e7a50.tar.bz2
LiterateLB-186bfae321f9d34b702a7c31b420eed4d10e7a50.tar.lz
LiterateLB-186bfae321f9d34b702a7c31b420eed4d10e7a50.tar.xz
LiterateLB-186bfae321f9d34b702a7c31b420eed4d10e7a50.tar.zst
LiterateLB-186bfae321f9d34b702a7c31b420eed4d10e7a50.zip
Small changes to channel example
-rw-r--r--lbm.org75
-rw-r--r--tangle/LLBM/sdf.h9
-rw-r--r--tangle/channel.cu (renamed from tangle/channel-with-sphere.cu)17
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>()),