|
| 1 | +/* The Computer Language Benchmarks Game |
| 2 | + https://salsa.debian.org/benchmarksgame-team/benchmarksgame/ |
| 3 | +
|
| 4 | + contributed by Brad Chamberlain |
| 5 | + derived from the Chapel version by Albert Sidelnik and myself |
| 6 | +
|
| 7 | + This version is the same as nbody3 but with a micro-optimization |
| 8 | + to 'advance' that removes an unnecessary **3. |
| 9 | +*/ |
| 10 | + |
| 11 | +use Math; // to get access to 'pi' |
| 12 | + |
| 13 | +config const n = 10000; // The number of timesteps to simulate |
| 14 | + |
| 15 | +param solarMass = 4 * pi * pi, |
| 16 | + daysPerYear = 365.24; |
| 17 | + |
| 18 | + |
| 19 | +// |
| 20 | +// a record for representing the position, velocity, and mass of |
| 21 | +// bodies in the solar system |
| 22 | +// |
| 23 | +record body { |
| 24 | + var pos: 3*real; |
| 25 | + var vel: 3*real; |
| 26 | + var mass: real; |
| 27 | +} |
| 28 | + |
| 29 | +// |
| 30 | +// the bodies that we'll be simulating |
| 31 | +// |
| 32 | +var bodies = (/* sun */ |
| 33 | + new body(mass = solarMass), |
| 34 | + |
| 35 | + /* jupiter */ |
| 36 | + new body(pos = ( 4.84143144246472090e+00, |
| 37 | + -1.16032004402742839e+00, |
| 38 | + -1.03622044471123109e-01), |
| 39 | + vel = ( 1.66007664274403694e-03 * daysPerYear, |
| 40 | + 7.69901118419740425e-03 * daysPerYear, |
| 41 | + -6.90460016972063023e-05 * daysPerYear), |
| 42 | + mass = 9.54791938424326609e-04 * solarMass), |
| 43 | + |
| 44 | + /* saturn */ |
| 45 | + new body(pos = ( 8.34336671824457987e+00, |
| 46 | + 4.12479856412430479e+00, |
| 47 | + -4.03523417114321381e-01), |
| 48 | + vel = (-2.76742510726862411e-03 * daysPerYear, |
| 49 | + 4.99852801234917238e-03 * daysPerYear, |
| 50 | + 2.30417297573763929e-05 * daysPerYear), |
| 51 | + mass = 2.85885980666130812e-04 * solarMass), |
| 52 | + |
| 53 | + /* uranus */ |
| 54 | + new body(pos = ( 1.28943695621391310e+01, |
| 55 | + -1.51111514016986312e+01, |
| 56 | + -2.23307578892655734e-01), |
| 57 | + vel = ( 2.96460137564761618e-03 * daysPerYear, |
| 58 | + 2.37847173959480950e-03 * daysPerYear, |
| 59 | + -2.96589568540237556e-05 * daysPerYear), |
| 60 | + mass = 4.36624404335156298e-05 * solarMass), |
| 61 | + |
| 62 | + /* neptune */ |
| 63 | + new body(pos = ( 1.53796971148509165e+01, |
| 64 | + -2.59193146099879641e+01, |
| 65 | + 1.79258772950371181e-01), |
| 66 | + vel = ( 2.68067772490389322e-03 * daysPerYear, |
| 67 | + 1.62824170038242295e-03 * daysPerYear, |
| 68 | + -9.51592254519715870e-05 * daysPerYear), |
| 69 | + mass = 5.15138902046611451e-05 * solarMass) |
| 70 | + ); |
| 71 | + |
| 72 | +param numBodies = bodies.size; // the number of bodies being simulated |
| 73 | + |
| 74 | +proc main() { |
| 75 | + initSun(); // initialize the sun's velocity |
| 76 | + |
| 77 | + writef("%.9r\n", energy()); // print the initial energy |
| 78 | + |
| 79 | + for 1..n do // simulate 'n' timesteps |
| 80 | + advance(0.01); |
| 81 | + |
| 82 | + writef("%.9r\n", energy()); // print the final energy |
| 83 | +} |
| 84 | + |
| 85 | +// |
| 86 | +// compute the sun's initial velocity |
| 87 | +// |
| 88 | +proc initSun() { |
| 89 | + const p = + reduce (for b in bodies do (b.vel * b.mass)); |
| 90 | + bodies[0].vel = -p / solarMass; |
| 91 | +} |
| 92 | + |
| 93 | +// |
| 94 | +// advance the positions and velocities of all the bodies |
| 95 | +// |
| 96 | +proc advance(dt) { |
| 97 | + for param i in 0..<numBodies { |
| 98 | + for param j in i+1..<numBodies { |
| 99 | + const dpos = bodies[i].pos - bodies[j].pos, |
| 100 | + dposNormSq = sumOfSquares(dpos), |
| 101 | + mag = dt / (dposNormSq * sqrt(dposNormSq)); |
| 102 | + |
| 103 | + bodies[i].vel -= dpos * bodies[j].mass * mag; |
| 104 | + bodies[j].vel += dpos * bodies[i].mass * mag; |
| 105 | + } |
| 106 | + } |
| 107 | + |
| 108 | + for param i in 0..<numBodies do |
| 109 | + bodies[i].pos += dt * bodies[i].vel; |
| 110 | +} |
| 111 | + |
| 112 | +// |
| 113 | +// compute the energy of the bodies |
| 114 | +// |
| 115 | +proc energy() { |
| 116 | + var e = 0.0; |
| 117 | + |
| 118 | + for i in 0..<numBodies { |
| 119 | + e += 0.5 * bodies[i].mass * sumOfSquares(bodies[i].vel); |
| 120 | + for j in i+1..<numBodies { |
| 121 | + e -= (bodies[i].mass * bodies[j].mass) |
| 122 | + / sqrt(sumOfSquares(bodies[i].pos - bodies[j].pos)); |
| 123 | + } |
| 124 | + } |
| 125 | + |
| 126 | + return e; |
| 127 | +} |
| 128 | + |
| 129 | +// |
| 130 | +// compute the sum of squares of a 3-tuple's elements |
| 131 | +// |
| 132 | +inline proc sumOfSquares((x,y,z)) { |
| 133 | + return x**2 + y**2 + z**2; |
| 134 | +} |
0 commit comments