Simulating 1-D Convection in Clojure — From Equations to Arrays
Earlier this year I gave a talk at the first online Scinoj Light Conference, sharing a ongoing project to port Computational Fluid Dynamics(CFD) learning materials from Python to Clojure.
In this post, I’ll demonstrate a simple one-dimensional linear convection simulation implemented in Clojure using Java primitive arrays. This example shows a glimpse into the porting effort, with a focus on expressing numerical simulations using only built-in Clojure functions.
The Equation
(This section won’t take up too much of your time…)
We’re going to simulate the 1-D linear convection equation:
\[\frac{\partial u }{\partial t} + c \frac{\partial u}{\partial x} = 0\]
This explains how the flow velocity u changes over time t and position x, with c which is the wave speed.
Instead of focusing on the physical interpretation(since I am no expert), the main focus on this post will be its implementation side - expressing it numerically and coding it in Clojure!
Using a finite-difference scheme, the discretized form becomes:
\[\frac{u_i^{n+1} - u_i^n}{\Delta t} + c \frac{u_i^n - u_{i-1}^n}{\Delta x} = 0\]
Then, solving for u_i^{n+1} gives:
\[u_i^{n+1} = n_i^n - c \frac{\Delta t}{\Delta x}(u_i^n - u_{i-1}^n)\]
Initial Setup
We begin by defining initial simulation parameters to start:
- nx: number of sliced steps for spatial point
xfromx-startandx-end - nt: number of time steps we want to propagate
- dx: each sliced
xcalculated fromdx = (x-end - x-start) / (nx - 1) - dt: sliced each time step
- c: speed of wave
(def init-params
{:x-start 0
:x-end 2
:nx 41
:nt 25
:dx (/ (- 2 0) (dec 41))
:dt 0.025
:c 1})Creating the x Grid
With the given init-params we defined earlier, we create a float-array of spatial points x:
(def array-x (let [{:keys [nx x-start x-end]} init-params
arr (float-array nx)
step (/ (- x-end x-start) (dec nx))]
(dotimes [i nx]
(aset arr i (float (* i step))))
arr))[0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6,
0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.05, 1.1, 1.15, 1.2,
1.25, 1.3, 1.35, 1.4, 1.45, 1.5, 1.55, 1.6, 1.65, 1.7, 1.75, 1.8,
1.85, 1.9, 1.95, 2.0]Defining Initial Flow Velocity Condition
The initial flow velocity(when t = 0) is 2 when x ∈ [0.5, 1.0], and is 1 elsewhere:
(def init-cond-fn #(float (if (and (>= % 0.5) (<= % 1.0)) 2 1)))(def array-u
(let [nx (:nx init-params)
u (float-array nx)]
(dotimes [i nx]
(let [x-i (aget array-x i)
u-i (init-cond-fn x-i)]
(aset u i u-i)))
u))(def init-array-u (float-array array-u))[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0,
2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]We can visualize this initial u:
Wait, Why dotimes?
Since we’re working with mutable Java arrays(float-array), dotimes is an efficient choice here. Because it gives direct, index-based iteration. And it pairs naturally with aget and aset for reading and writing array values.
Implementing and the Simulation
With the initial setup complete, we now apply the discretized convection equation at each time step.
Step Function
Given the previous time step’s array-u and the init-params, We compute and mutate the flow velocity in-place:
(defn mutate-linear-convection-u
[array-u {:keys [nx c dx dt]}]
(let [u_i (float-array array-u)]
(dotimes [i (- nx 2)]
(let [idx (inc i)
un-i (aget u_i idx)
un-i-1 (aget u_i i)
new-u-i (float (- un-i (* c (/ dt dx) (- un-i un-i-1))))]
(aset array-u idx new-u-i))))
array-u)Time Integration
We run the step function nt times to run our simulation over time.
(defn simulate!
[array-u {:keys [nt] :as init-params}]
(loop [n 0]
(if (= n nt)
array-u
(do (mutate-linear-convection-u array-u init-params) (recur (inc n))))))Finally, we visualize the resulting array-u:
(simulate! array-u init-params)[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0000007,
1.0000097, 1.0000782, 1.0004553, 1.0020387, 1.0073167, 1.0216427,
1.0538762, 1.1147615, 1.2121781, 1.3450189, 1.4999992, 1.6549712,
1.7877436, 1.8847833, 1.9440854, 1.9710407, 1.9710407, 1.9440854,
1.8847833, 1.7877436, 1.6549712, 1.4999992, 1.3450189, 1.2121781,
1.1147615, 1.0538762, 1.0216427, 1.0073167, 1.0]The plot shows how the flow velocity shifts from left to right over time, while also becoming smoother. Nice!
The Summary
This Simple example demonstrates a simulation process using low-level Java primitive arrays in Clojure.
Choosing this approach provided mutable, non-persistent data structures. While this deviates from idiomatic Clojure, it offers significant performance benefits for big-size numerical simulations. However, it comes with trade-offs; by opting for mutability, we give up the guarantees of immutability and structural sharing, making the code less safe and more error-prone.
As the porting project continues, we plan to evolve the design to better align with idiomatic Clojure principles. Stay tuned!
(comment
(require '[scicloj.clay.v2.api :as clay])
(clay/make! {:source-path "scicloj/cfd/intro/linear_1d_convection_with_array.clj"}))