examples: add an Lattice Boltzmann Method simulation (vlang#888)

Author: Bruno-Vdr

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
+}

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)))
+}

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.
+	}
+}

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]
+		}
+	}
+}

examples/lbm/profiles/asymetric.png

@@ -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.

examples/lbm/profiles/barrier.png

@@ -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: []
+}