Fortunately, Josh Barnes and Piet Hut had a better idea: recursively divide the simulation space into an oct-tree (oct- because we live in 3 dimensions, so splitting in two on each dimension gives 8 sub-trees for each tree) until each cell contains only one body. To compute the force on a particular body, iterate over the tree; for each cell, if the size of that cell is "small" compared with the distance to the body, approximate the forces due to all bodies in the cell by the force from a pseudo-body which resides at their center of mass and has a mass equal to the total mass of the bodies in the cell. (Whew---that's awkward to say. Hopefully it makes sense.) If the cell is not sufficiently "small", then consider its sub-cells recursively. In general (assuming you don't have to walk the entire tree to determine the force on each body), this process is O(log(n)) for each body's force, for a total time O(n*log(n)) to advance the system. Constructing the tree is also O(n*log(n)), so the whole process is O(n*log(n)). Much better!
Real cosmological codes are much more complicated than this prescription (in fact, they often use a different method for computing the really-long-range forces based on using Fourier transforms to solve the Poisson equation for the potential on a grid, but that's a story for another day). They also expend lots of work so they can use these trees in a distributed computation (if you want 1e9 particles, at 100 bytes/particle, just writing down your system state takes 100 GB, so you'd better distribute the computation or you'll not be able to fit it in main memory). I'm fortunately not involved in that type of work. When I say "simple" cosmological n-body code, I mean one which can simulate maybe 1e6 bodies on a workstation. The code I've posted below (released under the GPL) is a Barnes-Hut functor in OCaml. It forms the basis of my code, and, as you can see in the early comments, performs respectably, even with 1e6 bodies.
By the way, I'll leave it as an exercise for the reader to formulate the force accumulation algorithm in terms of fold_w_abort :). Also note that this code works (I think) in any number of dimensions (only tested in 3).
(** Barnes-Hut tree functor. Takes any structure which defines
[body_q] and [body_m] functions.
native code time for 1M-body tree construction: 82s. Memory usage
(just for construction): 223M. PowerBook G4 800 MHz, 1GB ram, Mac
OS 10.4.7, 18 August 2006.
Copyright (C) 2006 Will M. Farr
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
@author Will M. Farr *)
module type BODY =
sig
type body
val dim : int
val m : body -> float
val q : body -> float array
end
module Make(B : BODY) =
struct
(** [max a b] is a local, float-specialized, version of the
[Pervasives.max] function. *)
let max (a : float) (b : float) =
if a > b then a else b
(** [min a b] is a local, float-specialized, version of the
[Pervasives.min] function. *)
let min (a : float) (b : float) =
if a < b then a else b
(** Bounds are stored (for efficiency) as an array of [\[|low0, high0,
low1, high1, ... |\]]. Each bound is a half-open interval: [lowi]
<= xi < [highi]*)
type bounds = float array
(** A cell contains, in order, the total mass, center of mass, cell
bounds, cell size squared, sub-cells. *)
type tree =
| Empty
| Body of B.body
| Cell of float * float array * bounds * float * tree array
(** [m t] returns the mass of the tree [t] *)
let m = function
| Empty -> 0.0
| Body(b) -> B.m b
| Cell(m, _, _, _, _) -> m
(** [q t] returns the center of mass of the tree [t] *)
let q = function
| Empty -> Array.make 3 0.0
| Body(b) -> B.q b
| Cell(_, com, _, _, _) -> com
(** [low_bound bds i] returns the lower bound on dimension [i] given
bounds [bds]. *)
let low_bound (bds : bounds) i = bds.(2*i)
(** [high_bound bds i] return the high bound on dimension [i] given
bounds [bds]. *)
let high_bound (bds : bounds) i = bds.(2*i+1)
(** [in_bounds bs v] checks whether [v] is in the bounds given by
[bs]. *)
let in_bounds (bs : bounds) v =
let n = Array.length v in
let rec loop i =
if i >= n then
true
else if v.(i) >= low_bound bs i && v.(i) < high_bound bs i then
loop (i + 1)
else
false in
loop 0
(** [even i] and [odd i] *)
let even i = (i mod 2) = 0
let odd i = not (even i)
(** [pow_int i n] computes [i]^[n]. *)
let rec pow_int i n =
if n = 0 then
1
else if n = 1 then
i
else if even n then
let nn = n / 2 in
(pow_int i nn)*(pow_int i nn)
else
let nn = n / 2 in
i*(pow_int i nn)*(pow_int i nn)
(** [bit_set i n] is [true] if bit [i] of the integer [n] (bit 0 is
the least-significant bit of [n]) is 1 and [false] otherwise. *)
let bit_set i n =
n land (1 lsl i) > 0
(** [sub_bounds bs] returns an array of the 2^d (where d is the
dimension of the space---[bs] is 2*d elements long) sub-bounds of
the [bs] obtained by splitting the bounds in half on each
dimension. *)
let sub_bounds bs =
let nb = 2*B.dim and
nsb = 1 lsl B.dim in
Array.init nsb
(fun i ->
(* i is now encodes (bitwise) whether the bound in each
dimension is high or low *)
let sb = Array.make nb 0.0 in
for j = 0 to B.dim - 1 do
(* j labels the dimension *)
let mid = ((low_bound bs j) +. (high_bound bs j))/.2.0 in
if bit_set j i then
begin
sb.(2*j) <- mid;
sb.(2*j+1) <- high_bound bs j
end
else
begin
sb.(2*j) <- low_bound bs j;
sb.(2*j+1) <- mid
end
done;
sb)
(** [make_null_bounds ()] returns a fresh set of bounds which enclose
{b no} possible object. *)
let make_null_bounds () =
Array.init (2*B.dim)
(fun i ->
if even i then
infinity
else
neg_infinity)
(** [expand q] creates a value which is a bit bigger than [q] so that
[bounds_of_bodies bs] returns bounds which guarantee to enclose
[bs]. *)
let expand =
let factor = sqrt epsilon_float in
fun q ->
q +. (abs_float q)*.factor
(** [bounds_of_bodies bs] returns a bounds which completely enclose
the given bodies [bs]. *)
let bounds_of_bodies bs =
let bds = make_null_bounds () in
List.iter
(fun b ->
let bq = B.q b in
Array.iteri
(fun i q ->
bds.(2*i) <- min bds.(2*i) q;
bds.(2*i+1) <- max bds.(2*i+1) (expand q))
bq)
bs;
bds
(** [bounds_size_squared bds] returns the size squared of the given
bounds (i.e. the sum of squares of distances along each
dimension).*)
let bounds_size_squared bds =
let size = ref 0.0 and
n = (Array.length bds)/2 in
for i = 0 to n - 1 do
size := !size +. (Vector.square (bds.(2*i) -. bds.(2*i+1)))
done;
!size
(** [mass_and_com sts] returns the mass and center-of-mass of the [sts] *)
let mass_and_com sts =
let mass =
Array.fold_left
(fun mass t ->
mass +. (m t))
0.0
sts in
let com =
Array.fold_left
(fun com t ->
match t with
| Empty -> com
| _ ->
let mt = m t and
qt = q t in
for i = 0 to Array.length qt - 1 do
com.(i) <- com.(i) +. qt.(i)*.mt/.mass
done;
com)
(Array.make B.dim 0.0)
sts in
(mass, com)
(** [tree_of_body_list bs] constructs a tree which contains [bs]. *)
let rec tree_of_body_list = function
| [] -> Empty
| [b] -> Body(b)
| bs ->
let bds = bounds_of_bodies bs in
let s = bounds_size_squared bds in
let sub_bds = sub_bounds bds in
let sub_trees = Array.map (fun bd ->
let sub_bs = List.filter (fun b -> in_bounds bd (B.q b)) bs in
tree_of_body_list sub_bs)
sub_bds in
let m, com = mass_and_com sub_trees in
Cell(m, com, bds, s, sub_trees)
(** [tree_of_bodies bs] returns the tree which contains [bs]. *)
let tree_of_bodies bs =
let lbs = Array.to_list bs in
tree_of_body_list lbs
(** [fold fn start t] is the fundamental tree iterator. Alas, there
is no guarantee what the order of application of [fn] is. *)
let rec fold fn start t =
match t with
| Empty -> start
| Body(_) -> fn start t
| Cell(_, _, _, _, sts) ->
Array.fold_left
(fold fn) (fn start t) sts
(** [fold_w_abort fn start t] folds [fn] over [t] (with initial value
[start]). [fn] is applied to each tree-node before it is applied
to the sub-nodes. If [fn] returns [(false, value)] value is
returned as the result of the fold for this entire branch---the
recursion stops *)
let rec fold_w_abort fn start t =
match t with
| Empty -> start
| Body(_) -> snd (fn start t)
| Cell(_, _, _, _, sts) ->
let (cont, new_start) = fn start t in
if cont then
Array.fold_left (fold_w_abort fn) new_start sts
else
new_start
(** [contains b t] returns [true] if [t] contains [b]. *)
let rec contains b t =
let q = B.q b in
match t with
| Empty -> false
| Body(b2) -> b == b2
| Cell(_, _, bds, _, sts) ->
if in_bounds bds q then begin
let n = Array.length sts in
let rec loop i =
if i >= n then
false
else if contains b sts.(i) then
true
else
loop (i + 1) in
loop 0
end else
false
end
No comments:
Post a Comment