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/main.sml
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 193 - (view) (download)

1 : monnier 193 (* main.sml
2 :     *
3 :     * COPYRIGHT (c) 1993, AT&T Bell Laboratories.
4 :     *
5 :     * This is the main driver for the Barnse-HutN-body code.
6 :     *)
7 :    
8 :     functor Main (V : VECTOR) : sig
9 :    
10 :     structure S : SPACE
11 :     structure V : VECTOR
12 :     structure L : LOAD
13 :     structure G : GRAV
14 :    
15 :     val srand : int -> unit
16 :     (* reset the random number generator *)
17 :    
18 :     val utestdata : int -> S.body list
19 :     (* generate the uniform model data *)
20 :    
21 :     val testdata : int -> S.body list
22 :     (* generate the Plummer model data *)
23 :    
24 :     val go : {
25 :     output : {n2bcalc:int, nbccalc:int, nstep:int, selfint:int, tnow:real}
26 :     -> unit,
27 :     bodies : S.body list, tnow : real, tstop : real,
28 :     dtime : real, eps : real, tol : real,
29 :     rmin : real V.vec, rsize : real
30 :     } -> unit
31 :    
32 :     val doit : unit -> unit
33 :    
34 :     end = struct
35 :    
36 :     structure V = V
37 :     structure S = Space(V)
38 :     structure L = Load(S)
39 :     structure G = Grav(S)
40 :     structure DataIO = DataIO(S)
41 :    
42 :     (* some math utilities *)
43 :     val pi = 3.14159265358979323846
44 :     fun pow (x, 0.0) = 1.0
45 :     | pow (x, y) = exp (y * ln x)
46 :    
47 :     (* random numbers *)
48 :     local
49 :     val seed = ref 0.0
50 :     in
51 :     fun srand s = (seed := real s)
52 :     fun xrand (xl, xh) = let
53 :     val r = Random.random (! seed)
54 :     in
55 :     seed := r;
56 :     xl + (((xh - xl) * r) / 2147483647.0)
57 :     end
58 :     end (* local *)
59 :    
60 :     (* default parameter values *)
61 :     val defaults = [
62 :     (* file names for input/output *)
63 :     "in=", (* snapshot of initial conditions *)
64 :     "out=", (* stream of output snapshots *)
65 :    
66 :     (* params, used if no input specified, to make a Plummer Model*)
67 :     "nbody=128", (* number of particles to generate *)
68 :     "seed=123", (* random number generator seed *)
69 :    
70 :     (* params to control N-body integration *)
71 :     "dtime=0.025", (* integration time-step *)
72 :     "eps=0.05", (* usual potential softening *)
73 :     "tol=1.0", (* cell subdivision tolerence *)
74 :     "fcells=0.75", (* cell allocation parameter *)
75 :    
76 :     "tstop=2.0", (* time to stop integration *)
77 :     "dtout=0.25", (* data-output interval *)
78 :    
79 :     "debug=false", (* turn on debugging messages *)
80 :     "VERSION=1.0" (* JEB 06 March 1988 *)
81 :     ]
82 :    
83 :     (* pick a random point on a sphere of specified radius. *)
84 :     fun pickshell rad = let
85 :     fun pickvec () = let
86 :     val vec = V.tabulate (fn _ => xrand(~1.0, 1.0))
87 :     val rsq = V.dotvp(vec, vec)
88 :     in
89 :     if (rsq > 1.0)
90 :     then pickvec ()
91 :     else V.mulvs (vec, rad / sqrt(rsq))
92 :     end
93 :     in
94 :     pickvec ()
95 :     end
96 :    
97 :     fun utestdata n = let
98 :     val mfrac = 0.999 (* mass cut off at mfrac of total *)
99 :     val coeff = 4.0
100 :     val rn = real n
101 :     val rsc = (3.0 * pi) / 16.0
102 :     val vsc = sqrt(1.0 / rsc)
103 :     fun mkBodies (0, cmr, cmv, l) = let
104 :     (* offset bodies by normalized cm coordinates. Also, reverse
105 :     * the list to get the same order of bodies as in the C version.
106 :     *)
107 :     val cmr = V.divvs(cmr, rn)
108 :     val cmv = V.divvs(cmv, rn)
109 :     fun norm ([], l) = l
110 :     | norm ((p as S.Body{pos, vel, ...}) :: r, l) = (
111 :     pos := V.subv(!pos, cmr);
112 :     vel := V.subv(!vel, cmv);
113 :     norm (r, p::l))
114 :     in
115 :     norm (l, [])
116 :     end
117 :     | mkBodies (i, cmr, cmv, l) = let
118 :     val r = 1.0 / sqrt (pow(xrand(0.0, mfrac), ~2.0/3.0) - 1.0)
119 :     val pos = V.tabulate (fn _ => coeff*xrand(0.0, mfrac))
120 :     fun vN () = let (* von Neumann technique *)
121 :     val x = xrand(0.0,1.0)
122 :     val y = xrand(0.0,0.1)
123 :     in
124 :     if (y > x*x * (pow (1.0-x*x, 3.5))) then vN () else x
125 :     end
126 :     val v = ((sqrt 2.0) * vN()) / pow(1.0 + r*r, 0.25)
127 :     val vel = pickshell (vsc * v)
128 :     val body = S.Body{
129 :     proc = ref 0,
130 :     mass = 1.0 / rn,
131 :     pos = ref pos,
132 :     vel = ref vel,
133 :     acc = ref V.zerov,
134 :     phi = ref 0.0
135 :     }
136 :     in
137 :     mkBodies (i-1, V.addv(cmr, pos), V.addv(cmv, vel), body :: l)
138 :     end
139 :     in
140 :     mkBodies (n, V.zerov, V.zerov, [])
141 :     end (* testdata *)
142 :    
143 :     (* generate Plummer model initial conditions for test runs, scaled
144 :     * to units such that M = -4E = G = 1 (Henon, Hegge, etc).
145 :     * See Aarseth, SJ, Henon, M, & Wielen, R (1974) Astr & Ap, 37, 183.
146 :     *)
147 :     fun testdata n = let
148 :     val mfrac = 0.999 (* mass cut off at mfrac of total *)
149 :     val rn = real n
150 :     val rsc = (3.0 * pi) / 16.0
151 :     val vsc = sqrt(1.0 / rsc)
152 :     fun mkBodies (0, cmr, cmv, l) = let
153 :     (* offset bodies by normalized cm coordinates. Also, reverse
154 :     * the list to get the same order of bodies as in the C version.
155 :     *)
156 :     val cmr = V.divvs(cmr, rn)
157 :     val cmv = V.divvs(cmv, rn)
158 :     fun norm ([], l) = l
159 :     | norm ((p as S.Body{pos, vel, ...}) :: r, l) = (
160 :     pos := V.subv(!pos, cmr);
161 :     vel := V.subv(!vel, cmv);
162 :     norm (r, p::l))
163 :     in
164 :     norm (l, [])
165 :     end
166 :     | mkBodies (i, cmr, cmv, l) = let
167 :     val r = 1.0 / sqrt (pow(xrand(0.0, mfrac), ~2.0/3.0) - 1.0)
168 :     val pos = pickshell (rsc * r)
169 :     fun vN () = let (* von Neumann technique *)
170 :     val x = xrand(0.0,1.0)
171 :     val y = xrand(0.0,0.1)
172 :     in
173 :     if (y > x*x * (pow (1.0-x*x, 3.5))) then vN () else x
174 :     end
175 :     val v = ((sqrt 2.0) * vN()) / pow(1.0 + r*r, 0.25)
176 :     val vel = pickshell (vsc * v)
177 :     val body = S.Body{
178 :     proc = ref 0,
179 :     mass = 1.0 / rn,
180 :     pos = ref pos,
181 :     vel = ref vel,
182 :     acc = ref V.zerov,
183 :     phi = ref 0.0
184 :     }
185 :     in
186 :     mkBodies (i-1, V.addv(cmr, pos), V.addv(cmv, vel), body :: l)
187 :     end
188 :     in
189 :     mkBodies (n, V.zerov, V.zerov, [])
190 :     end (* testdata *)
191 :    
192 :     (* startup hierarchical N-body code. This either reads in or generates
193 :     * an initial set of bodies, and other parameters.
194 :     *)
195 :     fun startrun argv = let
196 :     val _ = GetParam.initParam(argv, defaults)
197 :     val {nbody, bodies, tnow, headline} = (case (GetParam.getParam "in")
198 :     of "" => let
199 :     val nbody = GetParam.getIParam "nbody"
200 :     in
201 :     if (nbody < 1)
202 :     then raise Fail "startrun: absurd nbody"
203 :     else ();
204 :     srand (GetParam.getIParam "seed");
205 :     { nbody = nbody,
206 :     bodies = testdata nbody,
207 :     tnow = 0.0,
208 :     headline = "Hack code: Plummer model"
209 :     }
210 :     end
211 :     | fname => DataIO.inputData fname
212 :     (* end case *))
213 :     in
214 :     { nbody = nbody,
215 :     bodies = bodies,
216 :     headline = headline,
217 :     outfile = GetParam.getParam "out",
218 :     dtime = GetParam.getRParam "dtime",
219 :     eps = GetParam.getRParam "eps",
220 :     tol = GetParam.getRParam "tol",
221 :     tnow = tnow,
222 :     tstop = GetParam.getRParam "tstop",
223 :     dtout = GetParam.getRParam "dtout",
224 :     debug = GetParam.getBParam "debug",
225 :     rmin = V.tabulate (fn _ => ~2.0),
226 :     rsize = 4.0
227 :     }
228 :     end
229 :    
230 :     (* advance N-body system one time-step. *)
231 :     fun stepSystem output {plist, dtime, eps, nstep, rmin, rsize, tnow, tol} = let
232 :     val dthf = 0.5 * dtime
233 :     val S.Space{rmin, rsize, root} = L.makeTree (plist, rmin, rsize)
234 :     (* recalculate accelaration *)
235 :     fun recalc ([], n2bcalc, nbccalc, selfint) = (n2bcalc, nbccalc, selfint)
236 :     | recalc (p::r, n2bcalc, nbccalc, selfint) = let
237 :     val S.Body{acc as ref acc1, vel, ...} = p
238 :     val {n2bterm, nbcterm, skipSelf} = G.hackGrav {
239 :     body = p, root = root, rsize = rsize, eps = eps, tol = tol
240 :     }
241 :     in
242 :     if (nstep > 0)
243 :     then (* use change in accel to make 2nd order *)
244 :     (* correction to vel. *)
245 :     vel := V.addv(!vel, V.mulvs(V.subv(!acc, acc1), dthf))
246 :     else ();
247 :     recalc (r, n2bcalc+n2bterm, nbccalc+nbcterm,
248 :     if skipSelf then selfint else (selfint+1))
249 :     end
250 :     (* advance bodies *)
251 :     fun advance (S.Body{pos, acc, vel, ...}) = let
252 :     val dvel = V.mulvs (!acc, dthf)
253 :     val vel1 = V.addv (!vel, dvel)
254 :     val dpos = V.mulvs (vel1, dtime)
255 :     in
256 :     pos := V.addv (!pos, dpos);
257 :     vel := V.addv (vel1, dvel)
258 :     end
259 :     val (n2bcalc, nbccalc, selfint) = recalc (plist, 0, 0, 0)
260 :     in
261 :     output {nstep=nstep, tnow=tnow, n2bcalc=n2bcalc, nbccalc=nbccalc, selfint=selfint};
262 :     app advance plist;
263 :     (nstep+1, tnow + dtime)
264 :     end
265 :    
266 :     (* given an initial configuration, run the simulation *)
267 :     fun go {
268 :     output, bodies, tnow, tstop,
269 :     dtime, eps, tol, rsize, rmin
270 :     } = let
271 :     val step = stepSystem output
272 :     fun loop (nstep, tnow) = if (tnow < tstop + (0.1 * dtime))
273 :     then loop (step {
274 :     plist = bodies, dtime = dtime, eps = eps, nstep = nstep,
275 :     rmin = rmin, rsize = rsize, tnow = tnow, tol = tol
276 :     })
277 :     else ()
278 :     in
279 :     loop (0, tnow)
280 :     end
281 :    
282 :     fun doit () = let
283 :     val { nbody, bodies, headline, outfile,
284 :     dtime, eps, tol, tnow, tstop, dtout,
285 :     debug, rsize, rmin
286 :     } = startrun []
287 :     fun output {nstep, tnow, n2bcalc, nbccalc, selfint} = DataIO.output{
288 :     bodies = bodies, nbody = nbody,
289 :     n2bcalc = n2bcalc, nbccalc = nbccalc,
290 :     selfint = selfint, tnow = tnow
291 :     }
292 :     in
293 :     DataIO.initOutput {
294 :     outfile = outfile, headline = headline, nbody = nbody, tnow = tnow,
295 :     dtime = dtime, eps = eps, tol = tol, dtout = dtout, tstop = tstop
296 :     };
297 :     go {
298 :     output=output, bodies=bodies, tnow=tnow, tstop=tstop,
299 :     dtime=dtime, eps=eps, tol=tol, rsize=rsize, rmin=rmin
300 :     };
301 :     DataIO.stopOutput()
302 :     end (* doit *)
303 :    
304 :     end; (* Main *)

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