The flow singularity approach is nice and easy but is sort of backwards. You specify the singularities and the resulting flow field is then computed. More often, you have some blob and want to know the resulting flow over it. To do this, a different approach is taken. You define the shape of your blob and other flow field parameters and then Laplace's equation is solved numerically. The solution is the resulting flow field, so this is called a "flow solver".
The trade off here is complexity. Not only for the underlying solution process (differential equation solver), but also for actually using the flow solver (grid generator, boundary conditions, etc.).
However, we found this really neat Python based solver:
This is super easy to use since the setup is simply a matrix. You fill the matrix with one of the possible flow elements, and then send it to the solver. The solver outputs the resulting flow field, also as a matrix.
This pairs nicely with the RGB matrix. The idea is to use the RGB matrix as a viewer for the flow solution. The solver won't actually run on the Matrix Portal* so we'll do something like this:
 Install the flow solver on a PC.
 Define our blob and run the solver on the PC.
 Output flow solution to a file.
 Copy solution file to Matrix Portal.
 Run viewer script on Matrix Portal.
*it probably could though, given enough time to port the code to CircuitPython.
Install Required Software
Do this software setup on the PC you want to run the flow solver on.
You'll need the following software. Click the links to go to installation instructions.
 Python  you may already have this installed.
 NumPy  used by the flow solver.
 Pillow/PIL  used to read BMP defining input geometry.
Then we need to install the flow solver software. If you are familiar with git, you can just clone the solver repo linked above. Otherwise, you can download the solver code as follows:

We'll work in the folder location where you saved the solver. So remember it.
Define Input Geometry
This is done with a BMP file. It should be a black and white image the same size as the RGB matrix. See the next section for details about creating this file.
Solver Script
Run this on your PC.
OK, now grab the solver script which will read in your BMP file, run the solver, and then output the results. To make things easy to run, save a copy of this file in the same folder as flow solver that was download above.
#====== # NOTE: Run this on your PC, not the Matrix Portal. #====== import sys import numpy as np from PIL import Image from ecoulements import systeme # load geometry grid = np.where(np.asarray(Image.open(sys.argv[1])), 1, 0) # add inlet / outlet flows inlet = np.array([2] * grid.shape[0]) outlet = np.array([3] * grid.shape[0]) grid = np.hstack((inlet[:, None], grid, outlet[:, None])) # add upper/ lower walls wall = np.array([0] * grid.shape[1]) grid = np.vstack((wall, grid, wall)) # solve _, VX, VY, _ = systeme.sol(grid) # save results to file OUTFILE = "flow_solution.py" with open(OUTFILE , "w") as fp: fp.write("nan = None\n") fp.write("solution = {\n") fp.write('"VX":\n') fp.write(str(VX[1:1, 1:1].tolist())) fp.write(',\n"VY":\n') fp.write(str(VY[1:1, 1:1].tolist())) fp.write("\n}\n") # done print("DONE! Results saved to", OUTFILE)
To run the solver, use the following command:
python flow_runner.py geometry.bmp
where geometry.bmp
is the BMP file you created to define your input geometry. You can use a different name if you want, but the BMP file should be located in the same folder location.
If it runs without errors, it should say DONE and the results will be in a new file named flow_solution.py
.
Move on to the next section for how to view those results!
Viewer Script
Run this on the Matrix Portal.
After running the solver script on your PC, you should end up with a file named flow_solution.py
. Copy that file to your CIRCUITPY folder. Then use the following script to view the results.
#====== # NOTE: Run this on the Matrix Portal. #====== import time import math import displayio from adafruit_matrixportal.matrix import Matrix from flow_solution import solution # User Config  SEEDS = ( # (x, y) starting location (0, 1), (0, 3), (0, 5), (0, 7), (0, 9), (0, 11), (0, 13), (0, 15), (0, 17), (0, 19), (0, 21), (0, 23), (0, 25), (0, 27), (0, 29), ) BACK_COLOR = 0x000000 # background fill SOLI_COLOR = 0xADAF00 # solids HEAD_COLOR = 0x00FFFF # leading particles TAIL_COLOR = 0x000A0A # trailing particles TAIL_LENGTH = 10 # length in pixels DELAY = 0.02 # smaller = faster # # use solution to define other items VX = solution['VX'] VY = solution['VY'] MATRIX_WIDTH = len(VX[0]) MATRIX_HEIGHT = len(VX) # meh...too lazy to list comp SOLIDS = [] for row in range(len(VX)): for col, v in enumerate(VX[row]): if v is None: SOLIDS.append((col, row)) # matrix and displayio setup matrix = Matrix(width=MATRIX_WIDTH, height=MATRIX_HEIGHT, bit_depth=6) display = matrix.display group = displayio.Group() display.show(group) bitmap = displayio.Bitmap(display.width, display.height, 4) palette = displayio.Palette(4) palette[0] = BACK_COLOR palette[1] = SOLI_COLOR palette[2] = HEAD_COLOR palette[3] = TAIL_COLOR tile_grid = displayio.TileGrid(bitmap, pixel_shader=palette) group.append(tile_grid) # global to store streamline data STREAMLINES = [] def compute_streamlines(): '''Compute streamline for each starting point (seed) defined.''' for seed in SEEDS: streamline = [] x, y = seed px = round(x) py = round(y) vx = VX[py][px] vy = VY[py][px] streamline.append( ((px, py), (vx, vy)) ) steps = 0 while x < MATRIX_WIDTH and steps < 2 * MATRIX_WIDTH: nx = round(x) ny = round(y) # if we've moved to a new pixel, store the info if nx != px or ny != py: streamline.append( ((nx, ny), (vx, vy)) ) px = nx py = ny if 0 <= nx < MATRIX_WIDTH and 0 <= ny < MATRIX_HEIGHT: vx = VX[ny][nx] vy = VY[ny][nx] if vx is None or vy is None: break x += vx y += vy steps += 1 # add streamline to global store STREAMLINES.append(streamline) def show_solids(): for s in SOLIDS: try: x, y = s bitmap[round(x), round(y)] = 1 except: # pylint: disable=bareexcept pass # just don't draw it def show_streamlines(): '''Draw the streamlines.''' for sl, head in enumerate(HEADS): try: streamline = STREAMLINES[sl] index = round(head) length = min(index, TAIL_LENGTH) # draw tail for data in streamline[indexlength:index]: x, y = data[0] bitmap[round(x), round(y)] = 3 # draw head bitmap[round(x), round(y)] = 2 except: # pylint: disable=bareexcept pass # just don't draw it def animate_streamlines(): '''Update the current location (head position) along each streamline.''' reset_heads = True for sl, head in enumerate(HEADS): # get associated streamline streamline = STREAMLINES[sl] # compute index index = round(head) # get velocity if index < len(streamline): vx, vy = streamline[index][1] reset_heads = False else: vx, vy = streamline[1][1] # move head HEADS[sl] += math.sqrt(vx*vx + vy*vy) if reset_heads: # all streamlines have reached the end, so reset to start for index, _ in enumerate(HEADS): HEADS[index] = 0 def update_display(): '''Update the matrix display.''' display.auto_refresh = False bitmap.fill(0) show_solids() show_streamlines() display.auto_refresh = True #========== # MAIN #========== print('Computing streamlines...', end='') compute_streamlines() print('DONE') HEADS = [0]*len(STREAMLINES) print('Flowing...') while True: animate_streamlines() update_display() time.sleep(DELAY)