Home My Page Projects Code Snippets Project Openings SML/NJ
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files

SCM Repository

[smlnj] Annotation of /sml/trunk/benchmarks/todo/barnes-but/grav.sml
ViewVC logotype

Annotation of /sml/trunk/benchmarks/todo/barnes-but/grav.sml

Parent Directory Parent Directory | Revision Log Revision Log


Revision 193 - (view) (download)

1 : monnier 193 (* grav.sml
2 :     *
3 :     * COPYRIGHT (c) 1993, AT&T Bell Laboratories.
4 :     *
5 :     * Gravity module for hierarchical N-body code; routines to compute gravity.
6 :     *)
7 :    
8 :     signature GRAV =
9 :     sig
10 :    
11 :     structure S : SPACE
12 :     structure V : VECTOR
13 :     sharing S.V = V
14 :    
15 :     val hackGrav : {body:S.body, root:S.node, rsize:real, tol:real, eps : real}
16 :     -> {n2bterm:int, nbcterm:int, skipSelf:bool}
17 :    
18 :     end; (* GRAV *)
19 :    
20 :     functor Grav (S : SPACE) : GRAV =
21 :     struct
22 :    
23 :     structure S = S
24 :     structure V = S.V
25 :    
26 :     fun walk {acc0, phi0, pos0, pskip, eps, rsize, tol, root} = let
27 :     val skipSelf = ref false
28 :     val nbcterm = ref 0 and n2bterm = ref 0
29 :     val tolsq = (tol * tol)
30 :     (* compute a single body-body or body-cell interaction. *)
31 :     fun gravsub (S.Empty, phi0, acc0, _) = (phi0, acc0)
32 :     | gravsub (p as S.Node{mass, pos, cell, ...}, phi0, acc0, memo) = let
33 :     val (dr, drsq) = (case memo
34 :     of NONE => let
35 :     val dr = V.subv(!pos, pos0)
36 :     in
37 :     (dr, V.dotvp(dr, dr) + (eps*eps))
38 :     end
39 :     | SOME(dr', drsq') => (dr', drsq' + (eps*eps))
40 :     (* end case *))
41 :     val phii = !mass / (sqrt drsq)
42 :     in
43 :     case cell
44 :     of (S.Cell _) => nbcterm := !nbcterm + 1
45 :     | _ => n2bterm := !n2bterm + 1
46 :     (* end case *);
47 :     (phi0 - phii, V.addv(acc0, V.mulvs(dr, phii / drsq)))
48 :     end (* gravsub *)
49 :     (* recursive routine to do hackwalk operation. This combines the
50 :     * subdivp and walksub routines from the C version.
51 :     *)
52 :     fun walksub (p, dsq, phi0, acc0) = (
53 :     (*print(implode[" walksub: dsq = ", makestring dsq, ", ", S.prNode p, "\n"]);*)
54 :     case p
55 :     of S.Empty => (phi0, acc0)
56 :     | (S.Node{cell = S.BodyCell body, ...}) =>
57 :     if (body = pskip)
58 :     then (skipSelf := true; (phi0, acc0))
59 :     else gravsub (p, phi0, acc0, NONE)
60 :     | (S.Node{cell = S.Cell a, pos, ...}) => let
61 :     val dr = V.subv(!pos, pos0)
62 :     val drsq = V.dotvp(dr, dr)
63 :     in
64 :     if ((tolsq * drsq) < dsq)
65 :     then let (* open p up *)
66 :     fun loop (i, phi0, acc0) = if (i < S.nsub)
67 :     then let
68 :     val (phi0', acc0') = walksub (
69 :     Array.sub(a, i), dsq/4.0, phi0, acc0)
70 :     in
71 :     loop (i+1, phi0', acc0')
72 :     end
73 :     else (phi0, acc0)
74 :     in
75 :     loop (0, phi0, acc0)
76 :     end
77 :     else gravsub (p, phi0, acc0, SOME(dr, drsq))
78 :     end
79 :     (* end case *))
80 :     val (phi0, acc0) = walksub (root, rsize*rsize, phi0, acc0)
81 :     in
82 :     { phi0 = phi0, acc0 = acc0,
83 :     nbcterm = !nbcterm, n2bterm = !n2bterm, skip = !skipSelf
84 :     }
85 :     end (* walk *)
86 :    
87 :     (* evaluate grav field at a given particle. *)
88 :     fun hackGrav {body as S.Body{pos, phi, acc, ...}, root, rsize, eps, tol} = let
89 :     val {phi0, acc0, nbcterm, n2bterm, skip} = walk {
90 :     acc0 = V.zerov, phi0 = 0.0, pos0 = !pos, pskip = body,
91 :     eps = eps, rsize = rsize, tol = tol, root = root
92 :     }
93 :     in
94 :     phi := phi0;
95 :     acc := acc0;
96 :     (**
97 :     app (fn (fmt, items) => print(Format.format fmt items)) [
98 :     ("pos = [%f %f %f]\n", map Format.REAL (V.explode(!pos))),
99 :     ("pos = [%f %f %f]\n", map Format.REAL (V.explode acc0)),
100 :     ("phi = %f\n", [Format.REAL phi0])
101 :     ];
102 :     raise Fail "";
103 :     **)
104 :     {nbcterm=nbcterm, n2bterm=n2bterm, skipSelf=skip}
105 :     end (* hackgrav *)
106 :    
107 :     end; (* Grav *)

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0