From fb8b76458f9f95597198b459484965ba4470c10c Mon Sep 17 00:00:00 2001
From: Bruno-Vdr <bruno@asterope.fr>
Date: Tue, 14 Jan 2025 12:12:23 +0100
Subject: [PATCH] examples: add an Lattice Boltzmann Method simulation (#888)
---
examples/lbm/cell.v | 176 +++++++++++++++
examples/lbm/colormap.v | 46 ++++
examples/lbm/lattice.v | 278 +++++++++++++++++++++++
examples/lbm/main.v | 330 ++++++++++++++++++++++++++++
examples/lbm/profiles/asymetric.png | Bin 0 -> 195 bytes
examples/lbm/profiles/barrier.png | Bin 0 -> 123 bytes
examples/lbm/profiles/circle.png | Bin 0 -> 249 bytes
examples/lbm/profiles/septum.png | Bin 0 -> 373 bytes
examples/lbm/profiles/test.png | Bin 0 -> 385 bytes
examples/lbm/readme.md | 32 +++
examples/lbm/v.mod | 7 +
11 files changed, 869 insertions(+)
create mode 100644 examples/lbm/cell.v
create mode 100644 examples/lbm/colormap.v
create mode 100644 examples/lbm/lattice.v
create mode 100644 examples/lbm/main.v
create mode 100644 examples/lbm/profiles/asymetric.png
create mode 100644 examples/lbm/profiles/barrier.png
create mode 100644 examples/lbm/profiles/circle.png
create mode 100644 examples/lbm/profiles/septum.png
create mode 100644 examples/lbm/profiles/test.png
create mode 100644 examples/lbm/readme.md
create mode 100644 examples/lbm/v.mod
diff --git a/examples/lbm/cell.v b/examples/lbm/cell.v
new file mode 100644
index 00000000..ded634e7
--- /dev/null
+++ b/examples/lbm/cell.v
@@ -0,0 +1,176 @@
+module main
+
+import math
+import rand
+
+// This file contains Code and data structures for a Lattice-Boltzmann-Method (LBM)
+// fluid flow simulation. The simulation is 2 Dimension, with 9 possible directions (D2Q9)
+//
+// 8 1 2
+//
+// 7 0 3
+//
+// 6 5 4
+
+// Vi is an enum for direction vector of D2Q9 Lattice.
+pub enum Vi {
+ center = 0
+ north = 1
+ north_east = 2
+ //
+ east = 3
+ south_east = 4
+ south = 5
+ //
+ south_west = 6
+ west = 7
+ north_west = 8
+}
+
+// vfmt off
+const opp = [
+ Vi.center, Vi.south, Vi.south_west,
+ Vi.west, Vi.north_west, Vi.north,
+ Vi.north_east, Vi.east, Vi.south_east,
+]!
+
+// Array defined here in order to loop over a Vi enum.
+// Warning: This array must be coherent with Vi enum order !
+const vi = [
+ Vi.center, Vi.north, Vi.north_east,
+ Vi.east, Vi.south_east, Vi.south,
+ Vi.south_west, Vi.west, Vi.north_west,
+]!
+
+// wi is the weight of mini-cells.
+// Warning: This array must be coherent with Vi enum order !
+const wi = [
+ f64(16.0 / 36.0), 4.0 / 36.0, 1.0 / 36.0,
+ 4.0 / 36.0, 1.0 / 36.0, 4.0 / 36.0,
+ 1.0 / 36.0, 4.0 / 36.0, 1.0 / 36.0,
+]!
+// vfmt on
+
+// opposite returns vector of opposite direction. Yes Enum can have methods in V.
+// Warning: This array must be coherent with Vi enum order !
+fn (v Vi) opposite() Vi {
+ return opp[int(v)]
+}
+
+// Discrete velocity vectors, by component, with respect to SDL display orientation (North is negative on y axis)
+// f64 and int are provided to avoid unnecessary casts.
+// Warning: These vectors must be coherent with Vi enum order !
+const dvx_f = [0.0, 0.0, 1.0, 1.0, 1.0, 0.0, -1.0, -1.0, -1.0]!
+const dvy_f = [0.0, -1.0, -1.0, 0.0, 1.0, 1.0, 1.0, 0.0, -1.0]!
+const dvx_i = [0, 0, 1, 1, 1, 0, -1, -1, -1]!
+const dvy_i = [0, -1, -1, 0, 1, 1, 1, 0, -1]!
+
+struct Cell {
+mut:
+ obstacle bool
+ sp [9]f64
+}
+
+// Cell.new built, and initializes a default Cell.
+// Warning: sp values must be coherent with Vi enum order !
+fn Cell.new(o bool) Cell {
+ return Cell{
+ obstacle: o
+ sp: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]!
+ } // ! means fixed size array ! Will change a day !
+}
+
+// Return a Cell with all mini-cell set to Zero
+fn Cell.zero(o bool) Cell {
+ return Cell{
+ obstacle: o
+ sp: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]!
+ }
+}
+
+// normalize perform a normalisation on the cell so that average desnsity is rho0
+fn (mut c Cell) normalize() {
+ rho := c.rho()
+ for mut v in c.sp {
+ v *= (rho0 / rho)
+ }
+}
+
+// randomize add some randomness on the cell.
+fn (mut c Cell) randomize(i f64) {
+ for v in vi {
+ r := rand.f64() * i
+ c.sum(v, r)
+ }
+}
+
+// get returns mini-cell depending on the given speed vector.
+fn (c &Cell) get(v Vi) f64 {
+ return c.sp[int(v)]
+}
+
+// set forces a mini-cell value given its speed vector.
+fn (mut c Cell) set(v Vi, value f64) {
+ c.sp[int(v)] = value
+}
+
+// increase (add value) to a mini-cell value given its speed vector.
+fn (mut c Cell) sum(v Vi, value f64) {
+ c.sp[int(v)] += value
+}
+
+// rho computes whole cell's density
+fn (c &Cell) rho() f64 {
+ mut sum := 0.0
+ for v in c.sp {
+ sum += v
+ }
+ assert math.is_nan(sum) == false
+ return sum
+}
+
+// ux computes x (horizontal) component of cell speed vector.
+fn (c &Cell) ux() f64 {
+ rho := c.rho()
+ r := 1.0 / rho * (c.sp[Vi.north_east] + c.sp[Vi.east] + c.sp[Vi.south_east] - c.sp[Vi.south_west] - c.sp[Vi.west] - c.sp[Vi.north_west])
+ assert math.is_nan(r) == false
+ return r
+}
+
+// uy computes y (vertical) component of cell speed vector.
+fn (c &Cell) uy() f64 {
+ rho := c.rho()
+ r := 1.0 / rho * (-c.sp[Vi.north] - c.sp[Vi.north_east] - c.sp[Vi.north_west] +
+ c.sp[Vi.south_east] + c.sp[Vi.south] + c.sp[Vi.south_west])
+ assert math.is_nan(r) == false
+ return r
+}
+
+// ux_no_rho computes x (horizontal) component of cell speed vector, when rho is already known and passed as param.
+fn (c &Cell) ux_no_rho(rho f64) f64 {
+ r := 1.0 / rho * (c.sp[Vi.north_east] + c.sp[Vi.east] + c.sp[Vi.south_east] - c.sp[Vi.south_west] - c.sp[Vi.west] - c.sp[Vi.north_west])
+ return r
+}
+
+// uy_no_rho computes y (vertical) component of cell speed vector, when rho is already known and passed as param.
+fn (c &Cell) uy_no_rho(rho f64) f64 {
+ r := 1.0 / rho * (-c.sp[Vi.north] - c.sp[Vi.north_east] - c.sp[Vi.north_west] +
+ c.sp[Vi.south_east] + c.sp[Vi.south] + c.sp[Vi.south_west])
+ return r
+}
+
+// equ computes result of equilibrium function.
+fn (c &Cell) equ(i Vi) f64 {
+ rho := c.rho()
+ ux := c.ux_no_rho(rho)
+ uy := c.uy_no_rho(rho)
+
+ t1 := 3.0 * (ux * dvx_f[i] + uy * dvy_f[i])
+ mut t2 := (ux * dvx_f[i] + uy * dvy_f[i])
+
+ t2 *= t2 // t2^2
+ t2 *= (9.0 / 2.0)
+ t3 := (3.0 / 2.0) * (ux * ux + uy * uy)
+ r := wi[i] * rho * (1.0 + t1 + t2 - t3)
+ return r
+}
diff --git a/examples/lbm/colormap.v b/examples/lbm/colormap.v
new file mode 100644
index 00000000..716dc457
--- /dev/null
+++ b/examples/lbm/colormap.v
@@ -0,0 +1,46 @@
+module main
+
+import arrays
+
+struct Colormap {}
+
+// build a linear color map ARGB8888 from color c1 to c2, of steps colors
+pub fn Colormap.build(c1 u32, c2 u32, steps int) []u32 {
+ assert steps > 1, 'Error, generating colormap needs at least two steps.'
+ mut cm := []u32{cap: steps}
+
+ // split c1 & c2 to RGB channels.
+ mut c1_r := f64((c1 & 0xFF0000) >> 16)
+ mut c1_g := f64((c1 & 0xFF00) >> 8)
+ mut c1_b := f64((c1 & 0xFF))
+
+ c2_r := f64((c2 & 0xFF0000) >> 16)
+ c2_g := f64((c2 & 0xFF00) >> 8)
+ c2_b := f64((c2 & 0xFF))
+
+ delta_r := (c2_r - c1_r) / f64(steps - 1)
+ delta_g := (c2_g - c1_g) / f64(steps - 1)
+ delta_b := (c2_b - c1_b) / f64(steps - 1)
+
+ cm << (0xFF000000 | c1) // Emit exact start color.
+ for _ in 1 .. steps - 1 {
+ c1_r += delta_r
+ c1_g += delta_g
+ c1_b += delta_b
+
+ // Recompose 3 channels to ARGB8888 color.
+ mut c := u32(0xFF_00_00_00)
+ c |= u32(c1_r) << 16
+ c |= u32(c1_g) << 8
+ c |= u32(c1_b)
+ cm << c
+ }
+ cm << (0xFF000000 | c2) // Emit exact end color.
+ return cm
+}
+
+// dual creates a color map with start color, intermediate color and final color, equaly sized.
+fn Colormap.dual(c1 u32, c2 u32, c3 u32, steps int) []u32 {
+ assert steps > 2, 'Error, generating dual-colormap needs at least three steps.'
+ return arrays.append(Colormap.build(c1, c2, steps / 2), Colormap.build(c2, c3, steps - (steps / 2)))
+}
diff --git a/examples/lbm/lattice.v b/examples/lbm/lattice.v
new file mode 100644
index 00000000..e03b8551
--- /dev/null
+++ b/examples/lbm/lattice.v
@@ -0,0 +1,278 @@
+module main
+
+import math
+
+pub struct Lattice {
+ w int
+ h int
+mut:
+ m []Cell
+}
+
+// Allocate x*y Cells Lattice
+pub fn Lattice.new(x int, y int, profile &u8) Lattice {
+ mut ret := Lattice{x, y, []Cell{cap: x * y}}
+
+ for i in 0 .. (x * y) {
+ e := unsafe { profile[i] != 0 }
+ mut c := Cell.new(e)
+ ret.m << c
+ }
+ return ret
+}
+
+pub fn (l Lattice) str() string {
+ return 'Lattice[${l.w}x${l.h}]'
+}
+
+// total_rho compute the total density on the Lattice. This value is conserved.
+pub fn (l Lattice) total_rho() f64 {
+ mut t := 0.0
+ for c in l.m {
+ if c.obstacle {
+ continue
+ }
+ t += c.rho()
+ }
+ return t
+}
+
+// clear the Lattice : File fields with zeros.
+pub fn (mut l Lattice) clear() {
+ unsafe { vmemset(l.m.data, 0, u32(l.m.len) * sizeof(Cell)) }
+}
+
+// add_flow create an artificial flow of i intensity in v direction
+// It impacts all lattice cells. It's usually a good thing to call normalize()
+// method after that.
+pub fn (mut l Lattice) add_flow(i f64, v Vi) {
+ for mut c in l.m {
+ if c.obstacle {
+ continue
+ }
+ c.sum(v, i)
+ }
+}
+
+// normalize normalizes all lattice cells against rho0 so that average density
+// is equal to rho0.
+pub fn (mut l Lattice) normalize() {
+ for mut c in l.m {
+ if c.obstacle {
+ continue
+ }
+ c.normalize()
+ }
+}
+
+// max_ux returns maximal horizontal speed in the whole lattice
+pub fn (l Lattice) max_ux() f64 {
+ mut ux := 0.0
+
+ for c in l.m {
+ if c.obstacle {
+ continue
+ }
+ u := c.ux()
+
+ if u > ux {
+ ux = u
+ }
+ }
+ return ux
+}
+
+// min_ux returns minimal horizontal speed in the whole lattice
+pub fn (l Lattice) min_ux() f64 {
+ mut ux := 0.0
+
+ for c in l.m {
+ if c.obstacle {
+ continue
+ }
+ u := c.ux()
+
+ if u < ux {
+ ux = u
+ }
+ }
+ return ux
+}
+
+// max_uy returns maximal horizontal speed in the whole lattice
+pub fn (l Lattice) max_uy() f64 {
+ mut uy := 0.0
+
+ for c in l.m {
+ if c.obstacle {
+ continue
+ }
+ u := c.uy()
+
+ if u > uy {
+ uy = u
+ }
+ }
+ return uy
+}
+
+// min_uy returns minimal horizontal speed in the whole lattice
+pub fn (l Lattice) min_uy() f64 {
+ mut uy := 0.0
+
+ for c in l.m {
+ if c.obstacle {
+ continue
+ }
+ u := c.uy()
+
+ if u < uy {
+ uy = u
+ }
+ }
+ return uy
+}
+
+// max_rho returns maximal cell density in the whole lattice
+pub fn (l Lattice) max_rho() f64 {
+ mut r := 0.0
+
+ for c in l.m {
+ if c.obstacle {
+ continue
+ }
+ rho := c.rho()
+
+ if rho > r {
+ r = rho
+ }
+ }
+ return r
+}
+
+// min_rho returns maximal cell density in the whole lattice
+pub fn (l Lattice) min_rho() f64 {
+ mut r := 1_000_000.0
+
+ for c in l.m {
+ if c.obstacle {
+ continue
+ }
+ rho := c.rho()
+
+ if rho < r {
+ r = rho
+ }
+ }
+ return r
+}
+
+// mean_rho returns the mean rho value over all lattice
+fn (l Lattice) mean_rho() f64 {
+ mut i := 0.0
+ for c in l.m {
+ if c.obstacle {
+ continue
+ }
+ i += c.rho()
+ }
+
+ return i / l.m.len
+}
+
+// vorticity computes vorticity at given position.
+fn (l Lattice) vorticity(x int, y int) f64 {
+ if x > 0 && x < l.w - 1 {
+ if y > 0 && y < l.h - 1 {
+ ind := y * l.w + x
+ omega := (l.m[ind + 1].uy() - l.m[ind - 1].uy()) - (l.m[ind + l.w].ux() - l.m[ind - l.w].ux())
+ return omega
+ }
+ }
+ return 0
+}
+
+// randomize add random noise everywhere given i intensity.
+// It's usually a good thing to call normalize() method after that.
+pub fn (mut l Lattice) randomize(i f64) {
+ for mut c in l.m {
+ if c.obstacle {
+ continue
+ }
+ c.randomize(i)
+ }
+}
+
+// move applies simulation movements handling borders and profile collisions.
+// Note that the rightmost column is handled differently, and outgoing mini-particle
+// are re-injected on the left.
+pub fn (l Lattice) move(mut output Lattice) {
+ mut index := 0
+ output.clear()
+
+ for y in 0 .. l.h {
+ for x in 0 .. l.w {
+ output.m[index].obstacle = l.m[index].obstacle // Copy src reachable state to output
+ for m in vi { // For this cell, for all direction vectors or mini particles
+ mini_part := l.m[index].get(m)
+ if dst_ind := l.reachable(x, y, m) {
+ output.m[dst_ind].sum(m, mini_part) // move mini-particle
+ } else {
+ output.m[index].sum(m.opposite(), mini_part) // rebound mini-particle, don't move but invert direction vector.
+ }
+ }
+ index++
+ }
+ }
+}
+
+// collide performs the most sensible step: Collision between mini-particles.
+pub fn (mut l Lattice) collide() {
+ for mut c in l.m {
+ if c.obstacle == false {
+ mut new_cell := Cell.zero(false)
+
+ for m in vi {
+ f := c.get(m)
+ feq := c.equ(m)
+ assert math.is_nan(feq) == false
+ assert math.is_nan(f) == false
+ new_cell.set(m, f - ((1.0 / tau) * (f - feq)))
+ }
+ c = new_cell
+ }
+ }
+}
+
+// reachable test if destination Cell (base on x,y and Direction vector)
+// is cross-able or reachable. Lattice edge is limiting, as the profile as
+// well, none is returned on these case. For reachable, index of the
+// reachable Cell (in Lattice) is returned, for an easy update.
+fn (l Lattice) reachable(x int, y int, v Vi) ?int {
+ assert x >= 0
+ assert y >= 0
+
+ // Add direction vector to current position to get destination position.
+ mut nx := x + dvx_i[int(v)]
+ ny := y + dvy_i[int(v)]
+
+ if ny < 0 || ny >= l.h {
+ return none
+ }
+
+ if nx < 0 {
+ nx = nx + l.w
+ }
+
+ if nx >= l.w {
+ nx = nx - l.w
+ }
+
+ ind := nx + (l.w * ny) // Get 1D index in lattice.
+
+ return if l.m[ind].obstacle {
+ none
+ } else {
+ ind // Destination cell is obstacle free. Return it's index.
+ }
+}
diff --git a/examples/lbm/main.v b/examples/lbm/main.v
new file mode 100644
index 00000000..7d8df0a8
--- /dev/null
+++ b/examples/lbm/main.v
@@ -0,0 +1,330 @@
+/**
+ This program is about D2Q9 Lattice Boltzmann Method:
+ See : https://en.wikipedia.org/wiki/Lattice_Boltzmann_methods
+
+ It's a pet project in order to use V language: https://vlang.io/
+
+ The simulation is single threaded, probably buggy and should not be used
+ for serious things. It's very sensible to tau parameter that should be
+ carefully set. This parameter is related to fluid viscosity, and is set so
+ that the fluid speed doesn't exceed a speed limit that breaks simulation.
+ Too narrow passage (speed increased) may reach this limit.
+
+ profiles files MUST be of the same size of defined width and weight and
+ should be 8bits per pixels. Every non zero value is considered as an
+ obstacle.
+
+ to compile the program from within source directory:
+
+ v -prod .
+
+ or if you want gcc as compiler:
+
+ v -prod -cc gcc .
+
+ SDL module must be installed: https://vpm.vlang.io/packages/sdl
+ and post install script executed, see link.
+
+ The simulation is quite slow, but would you like to slow it down, just
+ uncomment the sdl.delay(...) in file.
+
+
+
+ This program is released under MIT license.
+ */
+module main
+
+import sdl
+import sdl.image as img
+import os
+import time
+
+const tau = 0.6 // Relaxation time, related to fluid viscosity.
+const rho0 = 32.0 // Average normalised density.
+const width = 512
+const height = 128
+const len_in_bytes = width * height * sizeof(u32)
+const obstacle_color = u32(0xFF004000)
+const low_color = u32(0xFFFFFF00)
+const middle_color = u32(0xFF000000)
+const high_color = u32(0xFF00FFFF)
+
+// Type for rendering methods, used as parameter.
+type Renderer = fn (l Lattice, cm []u32, mut output []u32)
+
+const renderers = [vorticity, v_speed, h_speed, densities]
+const renderers_names = ['vorticity', 'vertical speed', 'horizontal speed', 'density']
+
+fn main() {
+ argv := os.args.len
+
+ if argv != 2 {
+ println('Usage: lbm profile_file.png')
+ println(' e.g: ./lbm profiles/circle.png')
+ println('During simulation press "v" to show different parameters.')
+ return
+ }
+
+ if sdl.init(sdl.init_video) < 0 {
+ eprintln('sdl.init() error: ${unsafe { cstring_to_vstring(sdl.get_error()) }}')
+ return
+ }
+
+ flags := u32(sdl.WindowFlags.opengl) | u32(sdl.WindowFlags.resizable)
+ window := sdl.create_window(c'Lattice Boltzmann Method [D2Q9]', sdl.windowpos_centered,
+ sdl.windowpos_centered, width * 2, height * 2, flags)
+
+ if window == sdl.null {
+ eprintln('sdl.create_window() error: ${unsafe { cstring_to_vstring(sdl.get_error()) }}')
+ return
+ }
+
+ r_flags := u32(sdl.RendererFlags.accelerated)
+ renderer := sdl.create_renderer(window, -1, r_flags)
+
+ if renderer == sdl.null {
+ eprintln('sdl.create_renderer() error: ${unsafe { cstring_to_vstring(sdl.get_error()) }}')
+ return
+ }
+
+ mut tex := sdl.create_texture(renderer, sdl.Format.argb8888, sdl.TextureAccess.streaming,
+ width, height)
+
+ if tex == sdl.null {
+ eprintln('sdl.create_texture() error: ${unsafe { cstring_to_vstring(sdl.get_error()) }}')
+ return
+ }
+
+ defer {
+ sdl.destroy_texture(tex)
+ sdl.destroy_renderer(renderer)
+ sdl.destroy_window(window)
+ }
+
+ profile := img.load(os.args[1].str)
+ if profile == sdl.null {
+ eprintln('Error trying to load profile .png file: ${unsafe { cstring_to_vstring(sdl.get_error()) }}')
+ return
+ }
+
+ // Check size compatibility.
+ if profile.w != width || profile.h != height {
+ eprintln('Error, "${os.args[1]}" profile image must match lbm lattice size : ${profile.w}x${profile.h}')
+ return
+ }
+
+ // Check profile is 1 byte / pixel.
+ if (profile.pitch / width) != 1 {
+ eprintln('Error profile file must be 1 byte per pixel')
+ return
+ }
+
+ // Build a colormap to be used
+ cm := Colormap.dual(low_color, middle_color, high_color, 384)
+
+ // Now create Lattices, with respect to loaded profile.
+ mut src := Lattice.new(width, height, profile.pixels)
+
+ src.add_flow(1.0, Vi.east)
+ src.randomize(0.2)
+ src.normalize()
+
+ mut dst := Lattice.new(width, height, profile.pixels)
+
+ // Allocate pixel buffer to draw in.
+ mut pixel_buffer := []u32{len: width * height} // Dyn array heap allocated.
+
+ mut should_close := false
+ mut frame := 0
+ mut render_method := vorticity
+
+ println('Showing vorticiy. Press "v" to show different parameters.')
+
+ for {
+ evt := sdl.Event{}
+ for 0 < sdl.poll_event(&evt) {
+ match evt.@type {
+ .quit {
+ should_close = true
+ }
+ .keydown {
+ key := unsafe { sdl.KeyCode(evt.key.keysym.sym) }
+ if key == sdl.KeyCode.escape {
+ should_close = true
+ break
+ }
+ // Show next view
+ if key == sdl.KeyCode.v {
+ mut i := renderers.index(render_method)
+ i++
+ if i >= renderers.len {
+ i = 0
+ }
+ render_method = renderers[i]
+ println('Rendering : ${renderers_names[i]}')
+ }
+ }
+ else {}
+ }
+ }
+ if should_close {
+ break
+ }
+
+ mut stop_watch := time.new_stopwatch()
+ src.move(mut dst)
+ dst.collide()
+ render_method(dst, cm, mut pixel_buffer) // render_method can point different method !
+ draw_colormap(cm, mut pixel_buffer)
+
+ // swap src and dst buffers.
+ tmp := src
+ src = dst
+ dst = tmp
+
+ blit_pixels(tex, pixel_buffer)
+ frame++
+ stop_watch.stop()
+ // println('Frame ${frame}, loop : ${stop_watch.elapsed().milliseconds()} milliseconds. ')
+
+ sdl.render_clear(renderer)
+ sdl.render_copy(renderer, tex, sdl.null, sdl.null)
+ sdl.render_present(renderer)
+ // sdl.delay(10)
+ }
+}
+
+// No bound checking here
+// blit_pixels from buffer to texture.
+@[direct_array_access]
+fn blit_pixels(t &sdl.Texture, data []u32) {
+ mut pixels := unsafe { voidptr(nil) }
+ mut pitch := int(0)
+
+ success := sdl.lock_texture(t, sdl.null, &pixels, &pitch)
+ if success < 0 {
+ panic('sdl.lock_texture error: ${sdl.get_error()}')
+ }
+ unsafe { vmemcpy(pixels, data.data, len_in_bytes) }
+ sdl.unlock_texture(t)
+}
+
+fn draw_colormap(cm []u32, mut data []u32) {
+ data[0] = 0xFF000000
+ data[width] = 0xFF000000
+ data[2 * width] = 0xFF000000
+ for i in 0 .. cm.len {
+ data[i + 1] = 0xFF000000
+ data[i + width + 1] = cm[i]
+ data[i + 1 + 2 * width] = 0xFF000000
+ }
+ data[cm.len] = 0xFF000000
+ data[width + cm.len] = 0xFF000000
+ data[2 * width + cm.len] = 0xFF000000
+}
+
+// densities is a Renderer type function
+fn densities(l Lattice, cm []u32, mut output []u32) {
+ mut ind := 0
+
+ min_rho := l.min_rho()
+ max_rho := l.max_rho()
+ linear := (max_rho - min_rho) / (cm.len - 1)
+
+ for c in l.m {
+ if c.obstacle == true {
+ output[ind] = obstacle_color
+ ind++
+ continue
+ }
+
+ rho := int((c.rho() - min_rho) / linear)
+ output[ind] = cm[rho]
+ ind++
+ }
+}
+
+// h_speed is a Renderer type function
+fn h_speed(l Lattice, cm []u32, mut output []u32) {
+ mut ind := 0
+
+ min_ux := l.min_ux()
+ max_ux := l.max_ux()
+ linear := (max_ux - min_ux) / (cm.len - 1)
+
+ for c in l.m {
+ if c.obstacle == true {
+ output[ind] = obstacle_color
+ ind++
+ continue
+ }
+
+ rho := int((c.ux() - min_ux) / linear)
+ output[ind] = cm[rho]
+ ind++
+ }
+}
+
+// h_speed is a Renderer type function
+fn v_speed(l Lattice, cm []u32, mut output []u32) {
+ mut ind := 0
+
+ min_uy := l.min_uy()
+ max_uy := l.max_uy()
+ linear := (max_uy - min_uy) / (cm.len - 1)
+
+ for c in l.m {
+ if c.obstacle == true {
+ output[ind] = obstacle_color
+ ind++
+ continue
+ }
+
+ rho := int((c.uy() - min_uy) / linear)
+ output[ind] = cm[rho]
+ ind++
+ }
+}
+
+// vorticity is a Renderer type function
+fn vorticity(l Lattice, cm []u32, mut output []u32) {
+ mut min := 0.0
+ mut max := 0.0
+ cm2 := u32(cm.len / 2)
+ mut vorticity_table := []f64{len: l.w * l.h}
+
+ for y in 0 .. l.h {
+ for x in 0 .. l.w {
+ out := (y * l.w) + x
+ if l.m[out].obstacle {
+ vorticity_table[out] = -1000000.0
+ } else {
+ v := l.vorticity(x, y)
+ vorticity_table[out] = v
+ if min > v {
+ min = v
+ }
+ if max < v {
+ max = v
+ }
+ }
+ }
+ }
+
+ linear := (max - min) / f64(cm.len - 1)
+
+ for ind, v in vorticity_table {
+ if v < -100.0 {
+ output[ind] = obstacle_color
+ } else {
+ mut id := cm2 + u32(v / linear)
+
+ if id < 0 {
+ id = 0
+ } else if id >= cm.len {
+ id = u32(cm.len - 1)
+ }
+ output[ind] = cm[id]
+ }
+ }
+}
diff --git a/examples/lbm/profiles/asymetric.png b/examples/lbm/profiles/asymetric.png
new file mode 100644
index 0000000000000000000000000000000000000000..cf6f6d3112e9dc10b1c9ef8f9de7476feec381fc
GIT binary patch
literal 195
zcmeAS@N?(olHy`uVBq!ia0y~yU;;838W@>@r29-~4<N-B;1l8sr2jK8NdNCu0<tPR
zT^vI)?!CRfk+(rX!X=T1b*ZwYRNpg=9*4}QH*d}I{-Icz->_bZeM$c8E0@m$bwEHz
zM@GWZ*#|z(aG1*JvAo(iR&2u3G)C7~-ORg7n05G%$g4=72-^Hzp=@GM=oc+<xz{Gn
gALad_nj175#6{~xc6`5E<O>q=boFyt=akR{08&;xiU0rr
literal 0
HcmV?d00001
diff --git a/examples/lbm/profiles/barrier.png b/examples/lbm/profiles/barrier.png
new file mode 100644
index 0000000000000000000000000000000000000000..6f7baf53226a402e538fef2596a726000d2a6a2c
GIT binary patch
literal 123
zcmeAS@N?(olHy`uVBq!ia0y~yU;;838W@>@r29-~4<N-B;1l8sr2jK8NdNCu0<sJ|
zT^vI)?!7(Z$jQLK!D8@MK1%tAFyA%?ZH5L7Tc$<_1py8gD3#8@uqWEq(~&W`cvc@!
OFN3G6pUXO@geCytB^yNm
literal 0
HcmV?d00001
diff --git a/examples/lbm/profiles/circle.png b/examples/lbm/profiles/circle.png
new file mode 100644
index 0000000000000000000000000000000000000000..ee52c9975fdc9ebfb45dd78cc5cb6bc93ec80ec2
GIT binary patch
literal 249
zcmeAS@N?(olHy`uVBq!ia0y~yU;;838W@>@r29-~4<N-B;1l8sr2jK8NdNCu0<zY5
zx;Tbp+<QCSkgwT+hoxCVqta>Kq05IZS8(|+Skh6D60206x$VXrqs#Y2ww1=+2Wnw>
zp#52sODXK-Uy(wOrs^ZXLOO!eBbe55ICGgf8Z1k=K5N0&M8AV8O9Nj&p44<&v%*E`
zY~ZF}EY89?vfRhiZv2llFg@7F;K*eqm%9JJ|E0Aovij4$=e?=g8}3w?vcxNM-JA%|
tt5c?GuigEwbyq};d2nnl$l4FWdl+U5X;&Uk74-t~Jzf1=);T3K0RYjKT6zEg
literal 0
HcmV?d00001
diff --git a/examples/lbm/profiles/septum.png b/examples/lbm/profiles/septum.png
new file mode 100644
index 0000000000000000000000000000000000000000..808e0b0aee0d1939296a00d6e80d79d3fb273c4b
GIT binary patch
literal 373
zcmeAS@N?(olHy`uVBq!ia0y~yU;;838W@>@r29-~4<N-B;1l8sr2jK8NdNCuVqjp@
z@^o<w$+-9S#>U(O3IeSUokgakFkkaHo-}dCMTw)zJ4zZWG}oGF{tyjpF+UI<@ttSR
z+&OdR{Fz-7@SxcGgKJ}=PGXRQK!=D&5nBph37ZI4NUX>Xb7o(D<|B=as~2#tKESw9
zFHpg|h^yq9D&GdPP%abC?`#S`qFa}mrPi<Ba7E2g;pBa-HCtAmy1#x0L(BZf^S2n(
z=BYF}EOKad5MkkBWodM2WNXymc=Fa^Nkig-ocYX$9hL%RFT9iI+CRmGD`K(3La<DL
zLf`>|0L=*+0(boU{wXywUDUjgxmS<tj=%{K4SrTOmRBtf!mK|MuWA;6oaANBDz}SQ
zszlmEQ=uw%Qt*QO`)qGx)jcmc&u3_6n*2`0h4G)*n>AdI4z#ZNa9X97S)SK7w!Ftj
Q3>bb4p00i_>zopr0PKp0_5c6?
literal 0
HcmV?d00001
diff --git a/examples/lbm/profiles/test.png b/examples/lbm/profiles/test.png
new file mode 100644
index 0000000000000000000000000000000000000000..4a8a74e778fef01b8e34d1401dae89b652452f3c
GIT binary patch
literal 385
zcmeAS@N?(olHy`uVBq!ia0y~yU;;838W@>@r29-~4<N-B;1l8sr2jK8NdNCuVqjo2
z^K@|x$+-9S#>QS2Lyp$O3CvDj3q-t<y>AC?V`&#()7E`%R>|M_Tck2OWNIhL-(bA-
zeAcJJy9aph9H?jXWL(}5m=NG#FhPWa%Yt<ZQ;>)U7Yl0-lNaMbhd>351Q7|504@(8
z>w$;?mj>$+4Uj?!Ru*R0#w88GK%JTj5*|REJxt9^O^r(&<}lu#@Q3Hg#_5lkR+$IS
zdHv_cZ=q*ZlIQoDohhArJ+7;`_vsz)kFE>{9rO<DuK$+yv6JC-or1$|R>n*0EU$PO
z>kSh`*&qFKh?8IOn)z?8gQfw{4Fy0aO=5CwY-<Qg2sof20Cb}Ws}z%0<Hyninh2gF
z<Kl*31)vKxxJ+2403GDU7|(F$_1=RL1`8Ay)E{8n)xiK!A<w-3X6TnWOI6JcfMLnt
M>FVdQ&MBb@0MVj?kpKVy
literal 0
HcmV?d00001
diff --git a/examples/lbm/readme.md b/examples/lbm/readme.md
new file mode 100644
index 00000000..ac95f716
--- /dev/null
+++ b/examples/lbm/readme.md
@@ -0,0 +1,32 @@
+This program is about D2Q9 Lattice Boltzmann Method:
+See : https://en.wikipedia.org/wiki/Lattice_Boltzmann_methods
+
+It's a pet project in order to use V language: https://vlang.io/
+
+The simulation is single threaded, probably buggy and should not be used
+for serious things. It's very sensible to tau parameter that should be
+carefully set. This parameter is related to fluid viscosity, and is set so
+that the fluid speed doesn't exceed a speed limit that breaks simulation.
+Too narrow passage (speed increased) may reach this limit.
+
+profiles files MUST be of the same size of defined width and weight and
+should be 8bits per pixels. Every non zero value is considered as an
+obstacle.
+
+to compile the program from within source directory:
+
+v -prod .
+
+or if you want gcc as compiler:
+
+v -prod -cc gcc .
+
+SDL module must be installed: https://vpm.vlang.io/packages/sdl
+and post install script executed, see link.
+
+The simulation is quite slow, but would you like to slow it down, just
+uncomment the sdl.delay(...) in file.
+
+
+
+This program is released under MIT license.
diff --git a/examples/lbm/v.mod b/examples/lbm/v.mod
new file mode 100644
index 00000000..5ed17dab
--- /dev/null
+++ b/examples/lbm/v.mod
@@ -0,0 +1,7 @@
+Module {
+ name: 'lbm'
+ description: 'Lattice Boltzmann simulation, used as toy projet to learn about V language'
+ version: '0.0.1'
+ license: 'MIT'
+ dependencies: []
+}