Journey to Vortex Shedding in the Browser: Part 2

How some object-oriented principles inspired me to fully refactor my codebase to fully remove all conditionals in my rendering loop and all the benefits that come with this new approach.

Context: I'm working on a live CFD simulation that runs in the browser, with a goal to display physically-realistic vortex shedding. See the first post for a full introduction.

Motivation for the Refactor

My code was quickly becoming a mess where each component was very tightly coupled to the implementation of different components. For example my equation assembly code had knowledge of the data structure used by my mesh, as shown in the following code where I'm comparing array indices with the structured mesh array's dimensions to check for boundaries.

// Some code from my original implementation
func (bc *BoundaryConditions) HandleBoundaryRHS(
		i, j, nx, ny int, aE, aW, aN, aS, dx, dy float64
	) float64 {
	rhs := 0.0

	if j == nx-1 { // East boundary
		rhs += bc.getRHSContribution(East, aE, dy)
	}
	if j == 0 { // West boundary
		rhs += bc.getRHSContribution(West, aW, dy)
	}
	if i == 0 { // North boundary
		rhs += bc.getRHSContribution(North, aN, dx)
	}
	if i == ny-1 { // South boundary
		rhs += bc.getRHSContribution(South, aS, dx)
	}

	return rhs
}

If I wanted to change to an unstructured mesh with a totally different data structure I'd have to change my entire codebase. Not only that but running the conditionals every loop to check for boundary conditions was a total waste - a cell's boundaries were not going to change between render loops and so that information should be precomputed. In short, my codebase looked a bit like this:

Crucially the issue with the way things interacted was that it wasn't through a well designed interface, but rather by messing with something else's internal data. For example I had code similar to the following.

// mesh.go
type Mesh struct {
	Phi   []float32
	NX NY int
	// etc...
}


// equations.go
type Equation struct {
	mesh *Mesh
	// etc...
}

func (eq *Equation) Assemble() {
	for index, cell := range eq.mesh.Phi {
		if index == eq.mesh.NX {
			// boundary stuff
		}
		// etc...
	}
}

Clearly the Assemble function is using the internal data of the meshing struct and so these two functionalities (meshing and equation assembly) become very tightly coupled, making change extremely difficult.

So how do I fix this mess?

I was hugely inspired by this excellent Sandi Metz talk , as I do a lot of OOP Ruby programming, and so I decided to translate the OOP principle of writing declarative code using polymorphism into Golang. If you're interested in learning more there's no way I can explain the theory better than Sandi in the above conference talk, so watch that first, but I'll give an example of how it changed my code below.

I think a lot of the power in OOP principles lies in the act of taking time to name objects and their behaviours properly, rather than sticking to a procedural framework for functionality. This simple exercise allows you to define clear boundaries and relationships and make your code far more declarative. The boundary condition logic, for example, is one such case where the procedure requires knowledge of different things: Is the cell in question touching a boundary? Which boundary? What boundary condition is set for that boundary? How should it effect the equations? Having to contend with all these is why the previous code was so tightly coupled - you do need information about the mesh, geometry (face areas and centroid distances), the field and the boundary conditions set by the user. My trick for disentangling this was to think:

Each Cell has a number of Neighbours that all provide some Contribution to the system. Therefore to assemble the equation each cell needs to get contributions from each of its neighbours.

This may seem like just semantics but notice there is no mention of boundary conditions or mesh details in this anthropomorphised account of the procedure. The cell doesn't need to know anything about whether its neighbours are boundaries or where they live, it just needs to know how to get the neighbours contributions. This allows us to write very declaritive code as in the following:

// field.go
func (f *Field) AssembleSystem(dt float32) *solver.System {
	nCells := len(f.cells)
	matrix := f.system.Matrix
	rhs := f.system.RHS

	for i, cell := range f.cells {
		diagonal := cell.timeCoeff / dt
		b := cell.timeCoeff / dt * cell.Phi

		// Notice there are no if/else statements
		for _, n := range cell.neighbours {
			diagonal += n.CouplingCoefficient()
			b += n.RHSContribution()
			n.ApplyOffDiagonalContribution(matrix, i, nCells)
		}
		matrix.Set(i, i, diagonal)
		rhs[i] = b
	}

	for i, src := range f.sources {
		matrix.Add(i, i, src.DiagonalContribution())
		rhs[i] += src.RHSContribution()
	}

	return f.system
}

How come we don't have to write any if statements above? Surely that logic (whether a neighbour is a boundary) needs to belong somewhere? The answer is being intentional about writing factory methods that figure that out and decide what sort of neighbours to pass to the cell. This approach, leveraging polymorphism with Go's interface system has changed the game for me, as I can separate the code into a "setup" phase, where factory functions populate my data structures with all the right sort of objects, and then my render loop (the "Hot loop") can just rely on a thing's behavioural contract, rather than having to figure out what sort of thing it is. This approach is called "Tell don't ask" and I've used to it fully eliminate any conditionals from my hot loop in my new refactor.

A quick explanation of polymorphism in Go

If you've not used Go's interfaces before or thought about duck typing, here's a quick example using my boundary condition code. We have a type of interface which says "any object that has these methods is also a type of this":

type Neighbour interface {
	CouplingCoefficient() float32
	ApplyOffDiagonalContribution(matrix *tensor.Matrix, cellIndex, nCells int)
	RHSContribution() float32
}

And then we have an internal neigbour type that implements Neighbour:

type InternalNeighbour struct {
	diffusionCoeff float32
	cellIndex      int
}

func (n *InternalNeighbour) CouplingCoefficient() float32 {
	return n.diffusionCoeff
}

func (n *InternalNeighbour) ApplyOffDiagonalContribution(matrix *tensor.Matrix, cellIndex, nCells int) {
	matrix.Set(cellIndex, n.cellIndex, -1.0*n.CouplingCoefficient())
}

func (n *InternalNeighbour) RHSContribution() float32 { return 0.0 }

But also a dirichlet neighour type that has the same methods:

type DirichletNeighbour struct {
	diffusionCoeff float32
	boundaryValue  float32
}

func (dn *DirichletNeighbour) CouplingCoefficient() float32 { return dn.diffusionCoeff }

func (dn *DirichletNeighbour) ApplyOffDiagonalContribution(matrix *tensor.Matrix, cellIndex, nCells int) {}

func (dn *DirichletNeighbour) RHSContribution() float32 {
	return dn.diffusionCoeff * dn.boundaryValue
}

This means that our Cell can have a strongly-typed field "neighbours" that is an array (well, technically a slice if you know Go) of type Neighbour and so all Cell "knows" about it's neighbours is that you can get contributions from them - all the logic is relegated to factory functions that figure out what type of neighbours each cell should have before the rendering loop starts. This is the core of polymorphism - an object is defined in terms of its functionality rather than its internal data or implementation details, allowing for very loose coupling of data structures.

Creating some general principles

Based on the success of the refactor for boundary conditoins,I've decided to adopt the following principles when writing future code for this project.

1. No condiitonals in my hot loop (the render loop)
2. Conditionals can only live in well defined factory functions
3. Components (renderer, mesh etc) must be as agnostic to other components as possible

So far these have been working great. Thinking deeply about the ways that demarcated blocks of code interact with each other has great benefit for adding more complex functionality cleanly and general quality of life for me as the developer. Also, being intentional about using factory functions makes it easy to reason about what belongs in the setup phase and the hot-loop, which gives free performance benefits like precomputing things that won't change to reduce the amount of float multiplication on every cycle.

The end result

Next steps

Before I use this new paradigm to add new features like convection terms or unstructured meshes, I really need to replace my Gauss-Seidel solver. This is because the tiny 21x11 grid above is running (in development) at < 10 FPS which is just abysmal. Having added a profiler (which you can check yourself in the website console, just press F12) I can see that my optimised equation assembly code is only taking an average of 50 nanoseconds per render. However, the system solving code is taking an average of ~150ms per render and so definitely needs to be optimised. I think my next job will be to learn about and implement a conjugate gradient method or something a lot more efficient.