diff options
| -rw-r--r-- | lbm.org | 75 | ||||
| -rw-r--r-- | tangle/LLBM/sdf.h | 9 | ||||
| -rw-r--r-- | tangle/channel.cu (renamed from tangle/channel-with-sphere.cu) | 17 |
3 files changed, 62 insertions, 39 deletions
@@ -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>()), |