Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

nbody demo: now with more working #292

Closed
wants to merge 3 commits into from
Closed

nbody demo: now with more working #292

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d693603
Updates to the nbody demo.
lkuper Mar 25, 2011
d920801
Merge remote branch 'graydon/master' into nbody-demo-2
lkuper Mar 25, 2011
ad1e44a
Got the nbody demo working, and hopefully faster, too!
lkuper Mar 26, 2011
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
225 changes: 120 additions & 105 deletions src/test/bench/shootout/nbody.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,13 +15,12 @@ fn main() {
let vec[Body.props] bodies = NBodySystem.MakeNBodySystem();

for (int n in inputs) {
// TODO: make #fmt handle floats?
log NBodySystem.energy(bodies);

let int i = 0;
while (i < n) {
bodies = NBodySystem.advance(bodies, 0.01);
i = i+1;
NBodySystem.advance(bodies, 0.01);
i += 1;
}
log NBodySystem.energy(bodies);
}
Expand All @@ -38,7 +37,6 @@ native "rust" mod rustrt {
mod NBodySystem {

fn MakeNBodySystem() -> vec[Body.props] {
// can't iterate over a record? how about a vector, then?
let vec[Body.props] bodies = vec(
// these each return a Body.props
Body.sun(),
Expand All @@ -51,41 +49,61 @@ mod NBodySystem {
let float py = 0.0;
let float pz = 0.0;

for (Body.props body in bodies) {
px += body.vx * body.mass;
py += body.vy * body.mass;
pz += body.vz * body.mass;
let int i = 0;
while (i < 5) {
px += bodies.(i).vx * bodies.(i).mass;
py += bodies.(i).vy * bodies.(i).mass;
pz += bodies.(i).vz * bodies.(i).mass;

i += 1;
}
bodies.(0) = Body.offsetMomentum(bodies.(0), px, py, pz);

// side-effecting
Body.offsetMomentum(bodies.(0), px, py, pz);

ret bodies;
}

fn advance(vec[Body.props] bodies, float dt) -> vec[Body.props] {
for (Body.props ibody in bodies) {
fn advance(vec[Body.props] bodies, float dt) -> () {

let Body.props iBody = ibody;

for (Body.props jbody in bodies) {
let float dx = iBody.x - jbody.x;
let float dy = iBody.y - jbody.y;
let float dz = iBody.z - jbody.z;
let int i = 0;
while (i < 5) {
let int j = i+1;
while (j < 5) {
let float dx = bodies.(i).x - bodies.(j).x;
let float dy = bodies.(i).y - bodies.(j).y;
let float dz = bodies.(i).z - bodies.(j).z;

let float dSquared = dx * dx + dy * dy + dz * dz;

let float distance;
rustrt.squareroot(dSquared, distance);
let float mag = dt / (dSquared * distance);

bodies.(i).vx -= dx * bodies.(j).mass * mag;
bodies.(i).vy -= dy * bodies.(j).mass * mag;
bodies.(i).vz -= dz * bodies.(j).mass * mag;

bodies.(j).vx += dx * bodies.(i).mass * mag;
bodies.(j).vy += dy * bodies.(i).mass * mag;
bodies.(j).vz += dz * bodies.(i).mass * mag;

j += 1;

}
}

for (Body.props body in bodies) {
body.x += dt * body.vx;
body.y += dt * body.vy;
body.z += dt * body.vz;
i += 1;
}

ret bodies;
i = 0;
while (i < 5) {

bodies.(i).x += dt * bodies.(i).vx;
bodies.(i).y += dt * bodies.(i).vy;
bodies.(i).z += dt * bodies.(i).vz;

i += 1;
}
}

fn energy(vec[Body.props] bodies) -> float {
Expand All @@ -95,116 +113,113 @@ mod NBodySystem {
let float distance;
let float e = 0.0;

for (Body.props ibody in bodies) {

// do we need this?
let Body.props iBody = ibody;

e += 0.5 * iBody.mass *
( iBody.vx * iBody.vx
+ iBody.vy * iBody.vy
+ iBody.vz * iBody.vz );

for (Body.props jbody in bodies) {

// do we need this?
let Body.props jBody = jbody;

dx = iBody.x - jBody.x;
dy = iBody.y - jBody.y;
dz = iBody.z - jBody.z;
let int i = 0;
while (i < 5) {
e += 0.5 * bodies.(i).mass *
( bodies.(i).vx * bodies.(i).vx
+ bodies.(i).vy * bodies.(i).vy
+ bodies.(i).vz * bodies.(i).vz );

let int j = i+1;
while (j < 5) {
dx = bodies.(i).x - bodies.(j).x;
dy = bodies.(i).y - bodies.(j).y;
dz = bodies.(i).z - bodies.(j).z;

rustrt.squareroot(dx*dx + dy*dy + dz*dz, distance);
e -= (iBody.mass * jBody.mass) / distance;
e -= (bodies.(i).mass * bodies.(j).mass) / distance;

j += 1;
}

i += 1;
}
ret e;
}

}
}

mod Body {
const float PI = 3.141592;
const float SOLAR_MASS = 39.478417; // was 4 * PI * PI originally

const float PI = 3.141592653589793;
const float SOLAR_MASS = 39.478417604357432; // was 4 * PI * PI originally
const float DAYS_PER_YEAR = 365.24;

type props = rec(float x,
float y,
float z,
float vx,
float vy,
float vz,
type props = rec(mutable float x,
mutable float y,
mutable float z,
mutable float vx,
mutable float vy,
mutable float vz,
float mass);

fn jupiter() -> Body.props {
// current limitation of the float lexer: decimal part has to
// fit into a 32-bit int.

let Body.props p;
p.x = 4.841431e+00;
p.y = -1.160320e+00;
p.z = -1.036220e-01;
p.vx = 1.660076e-03 * DAYS_PER_YEAR;
p.vy = 7.699011e-03 * DAYS_PER_YEAR;
p.vz = -6.904600e-05 * DAYS_PER_YEAR;
p.mass = 9.547919e-04 * SOLAR_MASS;
ret p;
ret rec(
mutable x = 4.84143144246472090e+00,
mutable y = -1.16032004402742839e+00,
mutable z = -1.03622044471123109e-01,
mutable vx = 1.66007664274403694e-03 * DAYS_PER_YEAR,
mutable vy = 7.69901118419740425e-03 * DAYS_PER_YEAR,
mutable vz = -6.90460016972063023e-05 * DAYS_PER_YEAR,
mass = 9.54791938424326609e-04 * SOLAR_MASS
);
}

fn saturn() -> Body.props {
let Body.props p;
p.x = 8.343366e+00;
p.y = 4.124798e+00;
p.z = -4.035234e-01;
p.vx = -2.767425e-03 * DAYS_PER_YEAR;
p.vy = 4.998528e-03 * DAYS_PER_YEAR;
p.vz = 2.304172e-05 * DAYS_PER_YEAR;
p.mass = 2.858859e-04 * SOLAR_MASS;
ret p;
ret rec(
mutable x = 8.34336671824457987e+00,
mutable y = 4.12479856412430479e+00,
mutable z = -4.03523417114321381e-01,
mutable vx = -2.76742510726862411e-03 * DAYS_PER_YEAR,
mutable vy = 4.99852801234917238e-03 * DAYS_PER_YEAR,
mutable vz = 2.30417297573763929e-05 * DAYS_PER_YEAR,
mass = 2.85885980666130812e-04 * SOLAR_MASS
);
}

fn uranus() -> Body.props {
let Body.props p;
p.x = 1.289436e+01;
p.y = -1.511115e+01;
p.z = -2.233075e-01;
p.vx = 2.964601e-03 * DAYS_PER_YEAR;
p.vy = 2.378471e-03 * DAYS_PER_YEAR;
p.vz = -2.965895e-05 * DAYS_PER_YEAR;
p.mass = 4.366244e-05 * SOLAR_MASS;
ret p;
ret rec(
mutable x = 1.28943695621391310e+01,
mutable y = -1.51111514016986312e+01,
mutable z = -2.23307578892655734e-01,
mutable vx = 2.96460137564761618e-03 * DAYS_PER_YEAR,
mutable vy = 2.37847173959480950e-03 * DAYS_PER_YEAR,
mutable vz = -2.96589568540237556e-05 * DAYS_PER_YEAR,
mass = 4.36624404335156298e-05 * SOLAR_MASS
);
}

fn neptune() -> Body.props {
let Body.props p;
p.x = 1.537969e+01;
p.y = -2.591931e+01;
p.z = 1.792587e-01;
p.vx = 2.680677e-03 * DAYS_PER_YEAR;
p.vy = 1.628241e-03 * DAYS_PER_YEAR;
p.vz = -9.515922e-05 * DAYS_PER_YEAR;
p.mass = 5.151389e-05 * SOLAR_MASS;
ret p;
ret rec(
mutable x = 1.53796971148509165e+01,
mutable y = -2.59193146099879641e+01,
mutable z = 1.79258772950371181e-01,
mutable vx = 2.68067772490389322e-03 * DAYS_PER_YEAR,
mutable vy = 1.62824170038242295e-03 * DAYS_PER_YEAR,
mutable vz = -9.51592254519715870e-05 * DAYS_PER_YEAR,
mass = 5.15138902046611451e-05 * SOLAR_MASS
);
}

fn sun() -> Body.props {
let Body.props p;
p.mass = SOLAR_MASS;
ret p;
ret rec(
mutable x = 0.0,
mutable y = 0.0,
mutable z = 0.0,
mutable vx = 0.0,
mutable vy = 0.0,
mutable vz = 0.0,
mass = SOLAR_MASS
);
}

fn offsetMomentum(Body.props props,
float px,
float py,
float pz) -> Body.props {

// TODO: should we create a new one or mutate the original?
let Body.props p = props;
p.vx = -px / SOLAR_MASS;
p.vy = -py / SOLAR_MASS;
p.vz = -pz / SOLAR_MASS;
ret p;
impure fn offsetMomentum(&Body.props props,
float px,
float py,
float pz) -> () {
props.vx = -px / SOLAR_MASS;
props.vy = -py / SOLAR_MASS;
props.vz = -pz / SOLAR_MASS;
}

}