From d78bc512d6f40bae64102a2f18fa085d3756eef8 Mon Sep 17 00:00:00 2001 From: John-John Tedro Date: Wed, 11 Oct 2023 04:16:46 +0200 Subject: [PATCH] Add nbodies benches --- benches/benches/nbodies.lua | 123 +++++++++++++++++++++++++++++ benches/benches/nbodies.rn | 150 ++++++++++++++++++++++++++++++++++++ 2 files changed, 273 insertions(+) create mode 100644 benches/benches/nbodies.lua create mode 100644 benches/benches/nbodies.rn diff --git a/benches/benches/nbodies.lua b/benches/benches/nbodies.lua new file mode 100644 index 000000000..9a07db9de --- /dev/null +++ b/benches/benches/nbodies.lua @@ -0,0 +1,123 @@ +-- The Computer Language Benchmarks Game +-- https://salsa.debian.org/benchmarksgame-team/benchmarksgame/ +-- contributed by Mike Pall +-- modified by Geoff Leyland + +local sqrt = math.sqrt + +local PI = 3.141592653589793 +local SOLAR_MASS = 4 * PI * PI +local DAYS_PER_YEAR = 365.24 +local bodies = { + { -- Sun + x = 0, + y = 0, + z = 0, + vx = 0, + vy = 0, + vz = 0, + mass = SOLAR_MASS + }, + { -- Jupiter + x = 4.84143144246472090e+00, + y = -1.16032004402742839e+00, + z = -1.03622044471123109e-01, + vx = 1.66007664274403694e-03 * DAYS_PER_YEAR, + vy = 7.69901118419740425e-03 * DAYS_PER_YEAR, + vz = -6.90460016972063023e-05 * DAYS_PER_YEAR, + mass = 9.54791938424326609e-04 * SOLAR_MASS + }, + { -- Saturn + x = 8.34336671824457987e+00, + y = 4.12479856412430479e+00, + z = -4.03523417114321381e-01, + vx = -2.76742510726862411e-03 * DAYS_PER_YEAR, + vy = 4.99852801234917238e-03 * DAYS_PER_YEAR, + vz = 2.30417297573763929e-05 * DAYS_PER_YEAR, + mass = 2.85885980666130812e-04 * SOLAR_MASS + }, + { -- Uranus + x = 1.28943695621391310e+01, + y = -1.51111514016986312e+01, + z = -2.23307578892655734e-01, + vx = 2.96460137564761618e-03 * DAYS_PER_YEAR, + vy = 2.37847173959480950e-03 * DAYS_PER_YEAR, + vz = -2.96589568540237556e-05 * DAYS_PER_YEAR, + mass = 4.36624404335156298e-05 * SOLAR_MASS + }, + { -- Neptune + x = 1.53796971148509165e+01, + y = -2.59193146099879641e+01, + z = 1.79258772950371181e-01, + vx = 2.68067772490389322e-03 * DAYS_PER_YEAR, + vy = 1.62824170038242295e-03 * DAYS_PER_YEAR, + vz = -9.51592254519715870e-05 * DAYS_PER_YEAR, + mass = 5.15138902046611451e-05 * SOLAR_MASS + } +} + +local function advance(bodies, nbody, dt) + for i=1,nbody do + local bi = bodies[i] + local bix, biy, biz, bimass = bi.x, bi.y, bi.z, bi.mass + local bivx, bivy, bivz = bi.vx, bi.vy, bi.vz + for j=i+1,nbody do + local bj = bodies[j] + local dx, dy, dz = bix-bj.x, biy-bj.y, biz-bj.z + local mag = sqrt(dx*dx + dy*dy + dz*dz) + mag = dt / (mag * mag * mag) + local bm = bj.mass*mag + bivx = bivx - (dx * bm) + bivy = bivy - (dy * bm) + bivz = bivz - (dz * bm) + bm = bimass*mag + bj.vx = bj.vx + (dx * bm) + bj.vy = bj.vy + (dy * bm) + bj.vz = bj.vz + (dz * bm) + end + bi.vx = bivx + bi.vy = bivy + bi.vz = bivz + bi.x = bix + dt * bivx + bi.y = biy + dt * bivy + bi.z = biz + dt * bivz + end +end + +local function energy(bodies, nbody) + local e = 0 + for i=1,nbody do + local bi = bodies[i] + local vx, vy, vz, bim = bi.vx, bi.vy, bi.vz, bi.mass + e = e + (0.5 * bim * (vx*vx + vy*vy + vz*vz)) + for j=i+1,nbody do + local bj = bodies[j] + local dx, dy, dz = bi.x-bj.x, bi.y-bj.y, bi.z-bj.z + local distance = sqrt(dx*dx + dy*dy + dz*dz) + e = e - ((bim * bj.mass) / distance) + end + end + return e +end + +local function offsetMomentum(b, nbody) + local px, py, pz = 0, 0, 0 + for i=1,nbody do + local bi = b[i] + local bim = bi.mass + px = px + (bi.vx * bim) + py = py + (bi.vy * bim) + pz = pz + (bi.vz * bim) + end + b[1].vx = -px / SOLAR_MASS + b[1].vy = -py / SOLAR_MASS + b[1].vz = -pz / SOLAR_MASS +end + +local N = 20000 +local nbody = #bodies + +offsetMomentum(bodies, nbody) +io.write( string.format("%0.9f",energy(bodies, nbody)), "\n") +for i=1,N do advance(bodies, nbody, 0.01) end +io.write( string.format("%0.9f",energy(bodies, nbody)), "\n") diff --git a/benches/benches/nbodies.rn b/benches/benches/nbodies.rn new file mode 100644 index 000000000..c91b3cbef --- /dev/null +++ b/benches/benches/nbodies.rn @@ -0,0 +1,150 @@ +const PI = 3.141592653589793; +const SOLAR_MASS = 4.0 * PI * PI; +const DAYS_PER_YEAR = 365.24; + +fn advance(bodies, dt) { + for i in 0..bodies.len() { + let bi = bodies[i]; + let bix = bi.x; + let biy = bi.y; + let biz = bi.z; + let bimass = bi.mass; + let bivx = bi.vx; + let bivy = bi.vy; + let bivz = bi.vz; + + for j in i + 1..bodies.len() { + let bj = bodies[j]; + let dx = bix - bj.x; + let dy = biy - bj.y; + let dz = biz - bj.z; + let mag = (dx * dx + dy * dy + dz * dz).sqrt(); + mag = dt / (mag * mag * mag); + let bm = bj.mass * mag; + bivx = bivx - (dx * bm); + bivy = bivy - (dy * bm); + bivz = bivz - (dz * bm); + bm = bimass * mag; + bj.vx = bj.vx + (dx * bm); + bj.vy = bj.vy + (dy * bm); + bj.vz = bj.vz + (dz * bm); + } + + bi.vx = bivx; + bi.vy = bivy; + bi.vz = bivz; + bi.x = bix + dt * bivx; + bi.y = biy + dt * bivy; + bi.z = biz + dt * bivz; + } +} + +fn energy(bodies) { + let e = 0.0; + + for i in 0..bodies.len() { + let bi = bodies[i]; + let vx = bi.vx; + let vy = bi.vy; + let vz = bi.vz; + let bim = bi.mass; + + e = e + (0.5 * bim * (vx * vx + vy * vy + vz * vz)); + + for j in i + 1..bodies.len() { + let bj = bodies[j]; + let dx = bi.x - bj.x; + let dy = bi.y - bj.y; + let dz = bi.z - bj.z; + let distance = (dx * dx + dy * dy + dz * dz).sqrt(); + e = e - ((bim * bj.mass) / distance); + } + } + + e +} + +fn offset_momentum(bodies) { + let px = 0.0; + let py = 0.0; + let pz = 0.0; + + for i in 0..bodies.len() { + let bi = bodies[i]; + let bim = bi.mass; + px = px + (bi.vx * bim); + py = py + (bi.vy * bim); + pz = pz + (bi.vz * bim); + } + + bodies[0].vx = -px / SOLAR_MASS; + bodies[0].vy = -py / SOLAR_MASS; + bodies[0].vz = -pz / SOLAR_MASS; +} + +pub fn main() { + let bodies = [ + // Sun + #{ + x: 0.0, + y: 0.0, + z: 0.0, + vx: 0.0, + vy: 0.0, + vz: 0.0, + mass: SOLAR_MASS, + }, + // Jupiter + #{ + x: 4.84143144246472090e+00, + y: -1.16032004402742839e+00, + z: -1.03622044471123109e-01, + vx: 1.66007664274403694e-03 * DAYS_PER_YEAR, + vy: 7.69901118419740425e-03 * DAYS_PER_YEAR, + vz: -6.90460016972063023e-05 * DAYS_PER_YEAR, + mass: 9.54791938424326609e-04 * SOLAR_MASS, + }, + // Saturn + #{ + x: 8.34336671824457987e+00, + y: 4.12479856412430479e+00, + z: -4.03523417114321381e-01, + vx: -2.76742510726862411e-03 * DAYS_PER_YEAR, + vy: 4.99852801234917238e-03 * DAYS_PER_YEAR, + vz: 2.30417297573763929e-05 * DAYS_PER_YEAR, + mass: 2.85885980666130812e-04 * SOLAR_MASS, + }, + // Uranus + #{ + x: 1.28943695621391310e+01, + y: -1.51111514016986312e+01, + z: -2.23307578892655734e-01, + vx: 2.96460137564761618e-03 * DAYS_PER_YEAR, + vy: 2.37847173959480950e-03 * DAYS_PER_YEAR, + vz: -2.96589568540237556e-05 * DAYS_PER_YEAR, + mass: 4.36624404335156298e-05 * SOLAR_MASS, + }, + // Neptune + #{ + x: 1.53796971148509165e+01, + y: -2.59193146099879641e+01, + z: 1.79258772950371181e-01, + vx: 2.68067772490389322e-03 * DAYS_PER_YEAR, + vy: 1.62824170038242295e-03 * DAYS_PER_YEAR, + vz: -9.51592254519715870e-05 * DAYS_PER_YEAR, + mass: 5.15138902046611451e-05 * SOLAR_MASS, + } + ]; + + let N = 20000; + let nbody = bodies.len(); + + offset_momentum(bodies); + println!("{}", energy(bodies)); + + for i in 1..N { + advance(bodies, 0.01); + } + + println!("{}", energy(bodies)); +}