Fortress Basic Simulation

Generic Include Statements

component particlesim
import List.{...}
export Executable

Objects and Overloaded Operators

object Vect3D(var x : RR64, var y : RR64, var z : RR64)
end

object Particle(var pos:Vect3D,var vel:Vect3D,var radius:RR64,var mass:RR64)
end

opr juxtaposition(v:Vect3D, n:Number): Vect3D = Vect3D(v.x n, v.y n, v.z n)

opr +(a:Vect3D,b:Vect3D): Vect3D = Vect3D(a.x + b.x, a.y + b.y, a.z + b.z)

opr DOT(a:Vect3D,b:Vect3D):RR64 = a.x b.x + a.y b.y + a.z b.z

magnitude(a:Vect3D):RR64 = SQRT (a.x a.x + a.y a.y + a.z a.z)

Collision Code

collisions(pop,size,timestep) = do
    var i : ZZ32 = 0
    while i < size do
        for j <- 0:(size-1) do
            var distVect: Vect3D = pop[j].pos + (pop[i].pos (-1.0))
            mag = magnitude(distVect) 

            if i NE j AND mag < (pop[i].radius + pop[j].radius) then
                atomic do
                    vcm = ((pop[i].vel pop[i].mass) + (pop[j].vel pop[j].mass)) (1 / (pop[i].mass + pop[j].mass))

                    pop[i].vel := pop[i].vel + (vcm (-1.0))
                    pop[j].vel := pop[j].vel + (vcm (-1.0))

                    distVect := distVect (1 / mag)

                    if (distVect DOT pop[i].vel GE 0) then
                        var vnorm: Vect3D = distVect (distVect DOT pop[i].vel)
                        var e : RR64 := magnitude(vnorm) 0.34
                        e := e^(-0.234)
                        if (e GT 1) then
                            e := 1
                        end

                        pop[i].vel := pop[i].vel + (vnorm ((-1.0) (1+e)))

                        vnorm := distVect (distVect DOT pop[j].vel)
                        pop[j].vel := pop[j].vel + (vnorm ((-1.0) (1+e)))
                    end

                    pop[i].vel := pop[i].vel + vcm
                    pop[j].vel := pop[j].vel + vcm
                end
            end
        end
        i += 1
    end
end

Gravity

This code was a little trickier to write. Because of the implicit parallelism of for loops, the outer iteration is implemented using a sequential while loop in order to ensure thread-safe performance. There isn't much else of note in this code except for the clear usage of adjacency for multiplication and extensive use of overloaded operators.

gravity(pop,size,timestep) = do
    var i : ZZ32 = 0
    while i < size do
        for j <- (i):(size-1) do
            if i NE j then
                var distVect: Vect3D = pop[i].pos + (pop[j].pos (-1.0))
                dist = magnitude(distVect) 
                softDist = dist + 0.0000001
                mag = timestep / (softDist softDist dist)

                pop[i].vel := pop[i].vel + (distVect (-1.0 pop[j].mass) mag)
                pop[j].vel := pop[j].vel + (distVect pop[i].mass mag)
            end
        end
        i += 1
    end
end

Stepper

This code is called every timestep in order to actually perform the simulation given the particle system.

step(pop,size,timestep) = do
    gravity(pop,size,timestep)
    collisions(pop,size,timestep)
    for i <- 0:(size-1) do
        var stepVector : Vect3D = pop[i].vel timestep
        pop[i].pos := pop[i].pos + stepVector
    end
end

The step function takes an array of particles that is referred to as a 'population', the total size of said array, and a timestep. It first performs all of the gravity forcing betweeen the particles, then handles all of the collisions for the current time step, and finally does the integrator 'drift' step. Note that the for structure is parallel everywhere in fortress, so the step of the simulation happens concurrently.

Main (a la C)

run() = do
    println "Starting particle simulation..."
    var i : ZZ32 = 0
    stepsize = 0.0001
    popsize = 2
    steps = 1000000
    pop : Particle[2] = [
          Particle(Vect3D(1.0,0.0,0.0) ,Vect3D(0.0,1.0,0.0),0.1,0.000001)
          Particle(Vect3D(0.0,0.0,0.0),Vect3D(0.0,0.0,0.0),0.5,1.0)]
    while i < steps do
        var j : ZZ32 = 0
        step(pop,popsize,stepsize)
        while j < popsize do
            println i " " j " " pop[j].pos.x " " pop[j].pos.y " " pop[j].pos.z " " pop[j].vel.x " " pop[j].vel.y " " pop[j].vel.z " " pop[j].radius
            j += 1
        end
        i += 1
    end
end

end particlesim