1 (* From the SML
/NJ benchmark suite
. *)
4 * COPYRIGHT (c
) 1993, AT
&T Bell Laboratories
.
6 * The abstract interface
of vectors
and matrices
in some dimension
.
13 val dim
: int (* dimension
of the vectors
*)
15 val tabulate
: (int -> 'a
) -> 'a vec
17 val equal
: real vec
* real vec
-> bool
19 val addv
: (real vec
* real vec
) -> real vec
20 val subv
: (real vec
* real vec
) -> real vec
21 val dotvp
: (real vec
* real vec
) -> real
22 val crossvp
: (real vec
* real vec
) -> real vec
23 val addvs
: (real vec
* real) -> real vec
24 val mulvs
: (real vec
* real) -> real vec
25 val divvs
: (real vec
* real) -> real vec
27 val mapv
: ('a
-> 'b
) -> 'a vec
-> 'b vec
28 val map3v
: (('a
* 'b
* 'c
) -> 'd
) -> ('a vec
* 'b vec
* 'c vec
) -> 'd vec
29 val foldv
: ('a
* 'b
-> 'b
) -> 'a vec
-> 'b
-> 'b
30 val format
: {lp
: string, sep
: string, rp
: string, cvt
: 'a
-> string}
32 val explode
: 'a vec
-> 'a list
33 val implode
: 'a list
-> 'a vec
35 type matrix (* matrices are always
real valued
*)
38 val addm
: (matrix
* matrix
) -> matrix
39 val outvp
: (real vec
* real vec
) -> matrix
44 * COPYRIGHT (c
) 1993, AT
&T Bell Laboratories
.
46 * The quad
/oct
-tree representation
of space
.
54 datatype body
= Body
of {
64 | Cell
of node Array
.array
74 datatype space
= Space
of {
80 val nsub
: int (* number
of sub cells
/ cell (2 ^ V
.dim
) *)
82 val putCell
: (cell
* int * node
) -> unit
83 val getCell
: (cell
* int) -> node
84 val mkCell
: unit
-> cell
85 val mkBodyNode
: body
-> node
86 val mkCellNode
: cell
-> node
87 val eqBody
: body
* body
-> bool
90 val dumpTree
: node
-> unit
91 val prBody
: body
-> string
92 val prNode
: node
-> string
96 functor Space (V
: VECTOR
) : SPACE
=
101 datatype body
= Body
of {
103 pos
: real V
.vec ref
,
104 vel
: real V
.vec ref
,
105 acc
: real V
.vec ref
,
111 | Cell
of node Array
.array
117 pos
: real V
.vec ref
,
121 datatype space
= Space
of {
127 fun eqBody(Body
{mass
,pos
,vel
,acc
,phi
},
128 Body
{mass
=m1
,pos
=p1
,vel
=v1
,acc
=a1
,phi
=h1
}) =
129 (Real.==(mass
, m1
) andalso Real.==(!phi
, !h1
)
130 andalso V
.equal(!pos
, !p1
) andalso V
.equal(!vel
, !v1
)
131 andalso V
.equal(!acc
, !a1
))
133 (* number
of sub cells per
cell (2 ^ V
.dim
) *)
134 val nsub
= Word.toInt(Word.<<(0w1
, Word.fromInt V
.dim
))
136 fun putCell (Cell a
, i
, nd
) = Array
.update(a
, i
, nd
)
137 fun getCell (Cell a
, i
) = Array
.sub(a
, i
)
138 fun mkCell () = Cell(Array
.array(nsub
, Empty
))
139 fun mkBodyNode (body
as Body
{pos
, mass
, ...}) = Node
{
140 cell
= BodyCell body
,
144 fun mkCellNode cell
= Node
{cell
= cell
, mass
= ref
0.0, pos
= ref V
.zerov
}
148 val rfmt
= Real.toString
149 val vfmt
= V
.format
{lp
="[", rp
="]", sep
=",", cvt
= rfmt
}
151 fun prBody (Body
{mass
, pos
, vel
, acc
, phi
}) = String.concat
[
156 ", phi=", rfmt(!phi
), "}"
158 fun prNode Empty
= "Empty"
159 |
prNode (Node
{mass
, pos
, cell
}) = let
160 val cell
= (case cell
161 of (Cell _
) => "Cell"
162 |
(BodyCell b
) => (*prBody b
*) "Body"
173 fun dumpTree tree
= let
174 fun printf items
= TextIO.output(TextIO.stdOut
, String.concat items
)
175 fun indent i
= StringCvt.padLeft #
" " (i
+1) ""
176 fun dump (node
, l
) = let
177 fun dump
' (Node
{cell
=Cell a
, ...}) = let
178 fun dump
'' i
= (dump(Array
.sub(a
, i
), l
+1); dump
''(i
+1))
180 (dump
'' 0) handle _
=> ()
185 StringCvt.padLeft #
" " 2 (Int.toString l
),
199 * COPYRIGHT (c
) 1993, AT
&T Bell Laboratories
.
201 * Code to build the tree from a list
of bodies
.
211 val makeTree
: (S
.body list
* real V
.vec
* real) -> S
.space
215 functor Load (S
: SPACE
) : LOAD
=
221 exception NotIntCoord
223 fun rshift (n
, k
) = Word.toInt(Word.~
>>(Word.fromInt n
, Word.fromInt k
))
225 val IMAX
= 0x20000000 (* 2^
29 *)
226 val IMAXrs1
= rshift(IMAX
, 1)
227 val rIMAX
= real IMAX
229 (* compute integerized coordinates
. Raises the NotIntCoord
exception,
230 * if rp is out
of bounds
.
232 fun intcoord (rp
, rmin
, rsize
) = let
233 val xsc
= V
.divvs (V
.subv(rp
, rmin
), rsize
)
234 fun cvt x
= if ((0.0 <= x
) andalso (x
< 1.0))
235 then floor(rIMAX
* x
)
236 else raise NotIntCoord
241 (* determine which subcell to select
. *)
242 fun subindex (iv
, l
) = let
243 fun aux (v
, (i
, k
)) = if (Word.andb(Word.fromInt v
, Word.fromInt l
) <> 0w0
)
244 then (i
+ rshift(S
.nsub
, k
+1), k
+1)
247 #
1 (V
.foldv aux
iv (0, 0))
250 (* enlarge cubical
"box", salvaging existing tree
structure. *)
251 fun expandBox (nd
as S
.Body
{pos
, ...}, box
as S
.Space
{rmin
, rsize
, root
}) = (
252 (intcoord (!pos
, rmin
, rsize
); box
)
253 handle NotIntCoord
=> let
254 val rmid
= V
.addvs (rmin
, 0.5 * rsize
)
255 val rmin
' = V
.map3v (fn (x
,y
,z
) =>
256 if x
< y
then z
- rsize
else z
) (!pos
, rmid
, rmin
)
257 val rsize
' = 2.0 * rsize
258 fun mksub (v
, r
) = let
259 val x
= intcoord (v
, rmin
', rsize
')
260 val k
= subindex (x
, IMAXrs1
)
261 val cell
= S
.mkCell ()
263 S
.putCell (cell
, k
, r
); cell
266 of S
.Empty
=> S
.Space
{rmin
=rmin
', rsize
=rsize
', root
=root
}
270 root
= S
.mkCellNode (mksub (rmid
, root
))
278 (* insert a single node into the tree
*)
279 fun loadTree (body
as S
.Body
{pos
=posp
, ...}, S
.Space
{rmin
, rsize
, root
}) = let
280 val xp
= intcoord (!posp
, rmin
, rsize
)
281 fun insert (S
.Empty
, _
) = S
.mkBodyNode body
282 |
insert (n
as S
.Node
{cell
=S
.BodyCell _
, pos
=posq
, ...}, l
) = let
283 val xq
= intcoord (!posq
, rmin
, rsize
)
284 val k
= subindex (xq
, l
)
288 insert (S
.mkCellNode a
, l
)
290 |
insert (n
as S
.Node
{cell
, ...}, l
) = let
291 val k
= subindex (xp
, l
)
292 val subtree
= insert (S
.getCell (cell
, k
), rshift(l
, 1))
294 S
.putCell (cell
, k
, subtree
);
298 S
.Space
{rmin
= rmin
, rsize
= rsize
, root
= insert (root
, IMAXrs1
)}
301 (* descend tree finding center
-of-mass coordinates
. *)
302 fun hackCofM S
.Empty
= ()
303 |
hackCofM (S
.Node
{cell
= S
.BodyCell _
, ...}) = ()
304 |
hackCofM (S
.Node
{cell
= S
.Cell subcells
, mass
, pos
}) = let
305 fun sumMass (i
, totMass
, cofm
) = if (i
< S
.nsub
)
306 then (case Array
.sub(subcells
, i
)
307 of S
.Empty
=> sumMass (i
+1, totMass
, cofm
)
308 |
(nd
as S
.Node
{mass
, pos
, ...}) => let
312 sumMass (i
+1, totMass
+ m
, V
.addv(cofm
, V
.mulvs(!pos
, m
)))
317 pos
:= V
.divvs(cofm
, totMass
))
319 sumMass (0, 0.0, V
.zerov
)
322 (* initialize tree
structure for hack force calculation
. *)
323 fun makeTree (bodies
, rmin
, rsize
) = let
324 fun build ([], space
) = space
325 |
build ((body
as S
.Body
{mass
, ...}) :: r
, space
) =
326 if Real.==(mass
, 0.0) then build (r
, space
)
328 val box
= expandBox (body
, space
)
329 val box
= loadTree(body
, box
)
332 val (space
as S
.Space
{root
, ...}) =
333 build (bodies
, S
.Space
{rmin
=rmin
, rsize
=rsize
, root
=S
.Empty
})
339 end; (* functor Load
*)
342 * COPYRIGHT (c
) 1993, AT
&T Bell Laboratories
.
344 * Gravity module for hierarchical N
-body code
; routines to compute gravity
.
354 val hackGrav
: {body
:S
.body
, root
:S
.node
, rsize
:real, tol
:real, eps
: real}
355 -> {n2bterm
:int, nbcterm
:int, skipSelf
:bool}
359 functor Grav (S
: SPACE
) : GRAV
=
365 fun walk
{acc0
, phi0
, pos0
, pskip
, eps
, rsize
, tol
, root
} = let
366 val skipSelf
= ref
false
367 val nbcterm
= ref
0 and n2bterm
= ref
0
368 val tolsq
= (tol
* tol
)
369 (* compute a single body
-body or body
-cell interaction
. *)
370 fun gravsub (S
.Empty
, phi0
, acc0
, _
) = (phi0
, acc0
)
371 |
gravsub (p
as S
.Node
{mass
, pos
, cell
, ...}, phi0
, acc0
, memo
) = let
372 val (dr
, drsq
) = (case memo
374 val dr
= V
.subv(!pos
, pos0
)
376 (dr
, V
.dotvp(dr
, dr
) + (eps
*eps
))
378 |
SOME(dr
', drsq
') => (dr
', drsq
' + (eps
*eps
))
380 val phii
= !mass
/ (Math
.sqrt drsq
)
383 of (S
.Cell _
) => nbcterm
:= !nbcterm
+ 1
384 | _
=> n2bterm
:= !n2bterm
+ 1
386 (phi0
- phii
, V
.addv(acc0
, V
.mulvs(dr
, phii
/ drsq
)))
388 (* recursive routine to
do hackwalk operation
. This combines the
389 * subdivp
and walksub routines from the C version
.
391 fun walksub (p
, dsq
, phi0
, acc0
) = (
392 (*print(implode
[" walksub: dsq = ", makestring dsq
, ", ", S
.prNode p
, "\n"]);*)
394 of S
.Empty
=> (phi0
, acc0
)
395 |
(S
.Node
{cell
= S
.BodyCell body
, ...}) =>
396 if S
.eqBody(body
, pskip
)
397 then (skipSelf
:= true; (phi0
, acc0
))
398 else gravsub (p
, phi0
, acc0
, NONE
)
399 |
(S
.Node
{cell
= S
.Cell a
, pos
, ...}) => let
400 val dr
= V
.subv(!pos
, pos0
)
401 val drsq
= V
.dotvp(dr
, dr
)
403 if ((tolsq
* drsq
) < dsq
)
404 then let (* open p up
*)
405 fun loop (i
, phi0
, acc0
) = if (i
< S
.nsub
)
407 val (phi0
', acc0
') = walksub (
408 Array
.sub(a
, i
), dsq
/4.0, phi0
, acc0
)
410 loop (i
+1, phi0
', acc0
')
416 else gravsub (p
, phi0
, acc0
, SOME(dr
, drsq
))
419 val (phi0
, acc0
) = walksub (root
, rsize
*rsize
, phi0
, acc0
)
421 { phi0
= phi0
, acc0
= acc0
,
422 nbcterm
= !nbcterm
, n2bterm
= !n2bterm
, skip
= !skipSelf
426 (* evaluate grav field at a given particle
. *)
427 fun hackGrav
{body
as S
.Body
{pos
, phi
, acc
, ...}, root
, rsize
, eps
, tol
} = let
428 val {phi0
, acc0
, nbcterm
, n2bterm
, skip
} = walk
{
429 acc0
= V
.zerov
, phi0
= 0.0, pos0
= !pos
, pskip
= body
,
430 eps
= eps
, rsize
= rsize
, tol
= tol
, root
= root
436 app (fn (fmt
, items
) => print(Format
.format fmt items
)) [
437 ("pos = [%f %f %f]\n", map Format
.REAL (V
.explode(!pos
))),
438 ("pos = [%f %f %f]\n", map Format
.REAL (V
.explode acc0
)),
439 ("phi = %f\n", [Format
.REAL phi0
])
443 {nbcterm
=nbcterm
, n2bterm
=n2bterm
, skipSelf
=skip
}
449 * COPYRIGHT (c
) 1993, AT
&T Bell Laboratories
.
451 * I
/O routines for export version
of hierarchical N
-body code
.
459 val inputData
: string -> {
461 bodies
: S
.body list
,
466 (* output routines
*)
468 outfile
: string, headline
: string, nbody
: int, tnow
: real,
469 dtime
: real, eps
: real, tol
: real, dtout
: real, tstop
: real
472 nbody
: int, bodies
: S
.body list
, n2bcalc
: int, nbccalc
: int,
473 selfint
: int, tnow
: real
475 val stopOutput
: unit
-> unit
479 functor DataIO (S
: SPACE
) : DATA_IO
=
482 structure SS
= Substring
486 val atoi
= valOf
o Int.scan
StringCvt.DEC SS
.getc
488 (* NOTE
: this really out to be implemented using the lazy IO streams
,
489 * but SML
/NJ doesn
't implement these correctly yet
.
491 fun inputData fname
= let
492 val strm
= TextIO.openIn fname
493 val buf
= ref(SS
.full
"")
494 fun getLn () = (case (TextIO.inputLine strm
)
495 of NONE
=> raise Fail
"inputData: EOF"
496 | SOME s
=> buf
:= SS
.full s
499 val buf
' = SS
.dropl
Char.isSpace (!buf
)
502 then (getLn(); skipWS())
506 val (n
, ss
) = atoi (skipWS ())
510 fun readReal () = let
511 val (r
, ss
) = valOf (Real.scan SS
.getc (skipWS()))
515 val nbody
= readInt()
516 val _
= if (nbody
< 1)
517 then raise Fail
"absurd nbody"
520 val _
= if (ndim
<> V
.dim
)
521 then raise Fail
"absurd ndim"
523 val tnow
= readReal()
526 |
loop (n
, l
) = loop (n
-1, f() :: l
)
530 fun readVec () = V
.implode (rev (iter readReal ndim
))
531 val massList
= iter readReal nbody
532 val posList
= iter readVec nbody
533 val velList
= iter readVec nbody
534 fun mkBodies ([], [], [], l
) = l
535 |
mkBodies (m
::r1
, p
::r2
, v
::r3
, l
) = let
544 mkBodies(r1
, r2
, r3
, b
::l
)
549 bodies
= mkBodies (massList
, posList
, velList
, []),
551 headline
= concat
["Hack code: input file ", fname
, "\n"]
556 val timer
= ref (Timer
.startCPUTimer ())
558 fun initTimer () = timer
:= Timer
.startCPUTimer()
560 val {usr
, sys
, ...} = Timer
.checkCPUTimer(!timer
)
563 (Time
.toReal totTim
) / 60.0
571 strm
: TextIO.outstream
573 val outState
= ref (NONE
: out_state option
)
575 fun fprintf (strm
, items
) = TextIO.output(strm
, String.concat items
)
576 fun printf items
= fprintf(TextIO.stdOut
, items
)
577 fun pad n s
= StringCvt.padLeft #
" " n s
578 fun fmtInt (wid
, i
) = pad
wid (Int.toString i
)
579 fun fmtReal (wid
, prec
, r
) = pad
wid (Real.fmt (StringCvt.FIX(SOME prec
)) r
)
580 fun fmtRealE (wid
, prec
, r
) = pad
wid (Real.fmt (StringCvt.SCI(SOME prec
)) r
)
582 fun itemFmt r
= fmtReal (9, 4, r
)
583 val fmt
= V
.format
{lp
="", sep
="", rp
="", cvt
=itemFmt
}
585 fun printvec (init
, vec
) = printf
[
586 "\t ", pad
9 init
, fmt vec
, "\n"
590 fun stopOutput () = (case (! outState
)
592 |
(SOME
{strm
, ...}) => (TextIO.closeOut strm
; outState
:= NONE
)
595 fun initOutput
{outfile
, headline
, nbody
, tnow
, dtime
, eps
, tol
, dtout
, tstop
} = (
597 printf
["\n\t\t", headline
, "\n\n"];
598 printf (map (pad
12) ["nbody", "dtime", "eps", "tol", "dtout", "tstop"]);
600 printf
[fmtInt(12, nbody
), fmtReal(12, 5, dtime
)];
602 fmtInt(12, nbody
), fmtReal(12, 5, dtime
),
603 fmtReal(12, 4, eps
), fmtReal(12, 2, tol
),
604 fmtReal(12, 3, dtout
), fmtReal(12, 2, tstop
), "\n\n"
607 of "" => stopOutput()
608 | _
=> outState
:= SOME
{
612 strm
= TextIO.openOut outfile
616 (* compute set
of dynamical diagnostics
. *)
617 fun diagnostics bodies
= let
618 fun loop ([], arg
) = {
619 mtot
= #totM arg
, (* total mass
*)
620 totKE
= #totKE arg
, (* total kinetic energy
*)
621 totPE
= #totPE arg
, (* total potential energy
*)
622 cOfMPos
= #cOfMPos arg
, (* center
of mass
: position
*)
623 cOfMVel
= #cOfMVel arg
, (* center
of mass
: velocity
*)
624 amVec
= #amVec
arg (* angular momentum vector
*)
627 mass
, pos
=ref pos
, vel
=ref vel
, acc
=ref acc
, phi
=ref phi
629 val velsq
= V
.dotvp(vel
, vel
)
630 val halfMass
= 0.5 * mass
631 val posXmass
= V
.mulvs(pos
, mass
)
634 totM
= (#totM arg
) + mass
,
635 totKE
= (#totKE arg
) + halfMass
* velsq
,
636 totPE
= (#totPE arg
) + halfMass
* phi
,
637 keTen
= V
.addm(#keTen arg
, V
.outvp(V
.mulvs(vel
, halfMass
), vel
)),
638 peTen
= V
.addm(#peTen arg
, V
.outvp(posXmass
, acc
)),
639 cOfMPos
= V
.addv(#cOfMPos arg
, posXmass
),
640 cOfMVel
= V
.addv(#cOfMVel arg
, V
.mulvs(vel
, mass
)),
641 amVec
= V
.addv(#amVec arg
, V
.mulvs(V
.crossvp(pos
, vel
), mass
))
646 totM
= 0.0, totKE
= 0.0, totPE
= 0.0,
647 keTen
= V
.zerom
, peTen
= V
.zerom
,
648 cOfMPos
= V
.zerov
, cOfMVel
= V
.zerov
,
651 end (* diagnostics
*)
653 fun outputData (strm
, tnow
, nbody
, bodies
) = let
654 fun outInt i
= fprintf(strm
, [" ", Int.toString i
, "\n"])
655 fun outReal r
= fprintf(strm
, [" ", fmtRealE(21, 14, r
), "\n"])
656 fun prReal r
= fprintf(strm
, [" ", fmtRealE(21, 14, r
)])
658 fun out
[] = TextIO.output(strm
, "\n")
659 |
out (x
::r
) = (prReal x
; out r
)
667 app (fn (S
.Body
{mass
, ...}) => outReal mass
) bodies
;
668 app (fn (S
.Body
{pos
, ...}) => outVec(!pos
)) bodies
;
669 app (fn (S
.Body
{vel
, ...}) => outVec(!vel
)) bodies
;
670 printf
["\n\tparticle data written\n"]
673 fun output
{nbody
, bodies
, n2bcalc
, nbccalc
, selfint
, tnow
} = let
674 val nttot
= n2bcalc
+ nbccalc
675 val nbavg
= floor(real n2bcalc
/ real nbody
)
676 val ncavg
= floor(real nbccalc
/ real nbody
)
677 val data
= diagnostics bodies
680 printf (map (pad
9) [
681 "tnow", "T+U", "T/U", "nttot", "nbavg", "ncavg", "selfint",
686 fmtReal(9, 3, tnow
), fmtReal(9, 4, #totKE data
+ #totPE data
),
687 fmtReal(9, 4, #totKE data
/ #totPE data
), fmtInt(9, nttot
),
688 fmtInt(9, nbavg
), fmtInt(9, ncavg
), fmtInt(9, selfint
),
689 fmtReal(9, 2, cputime()), "\n\n"
691 printvec ("cm pos", #cOfMPos data
);
692 printvec ("cm vel", #cOfMVel data
);
693 printvec ("am pos", #amVec data
);
696 |
(SOME
{tout
, dtout
, dtime
, strm
}) =>
697 if ((tout
- 0.01 * dtime
) <= tnow
)
699 outputData (strm
, tnow
, nbody
, bodies
);
701 tout
=tout
+dtout
, dtout
=dtout
, dtime
=dtime
, strm
=strm
711 * COPYRIGHT (c
) 1993, AT
&T Bell Laboratories
.
714 structure GetParam
: sig
718 val initParam
: (string list
* string list
) -> unit
719 val getParam
: string -> string
720 val getIParam
: string -> int
721 val getRParam
: string -> real
722 val getBParam
: string -> bool
728 val defaults
= ref ([] : string list
)
730 (* ignore arg vector
, remember defaults
. *)
731 fun initParam (argv
, defl
) = defaults
:= defl
734 TextIO.output(TextIO.stdOut
, String.concat items
);
735 TextIO.flushOut
TextIO.stdOut
)
737 structure SS
= Substring
739 (* export version prompts user for value
. *)
740 fun getParam name
= let
741 fun scanBind
[] = NONE
742 |
scanBind (s
::r
) = let
743 val (_
, suffix
) = SS
.position
name (SS
.full s
)
745 if (SS
.isEmpty suffix
)
747 else SOME(SS
.string(SS
.triml (size name
+1) suffix
))
749 fun get default
= (case (TextIO.inputLine
TextIO.stdIn
)
751 | SOME
"\n" => default
752 | SOME s
=> substring(s
, 0, size s
- 1)
755 if (null (! defaults
))
756 then raise Fail
"getParam called before initParam"
758 case (scanBind (! defaults
))
760 prompt
["enter ", name
, " [", s
, "]: "];
762 | NONE
=> (prompt
["enter ", name
, ": "]; get
"")
769 fun get () = (case getParam name
of "" => get () | s
=> s
)
772 (valOf (scanFn param
)) handle _
=> (cvt
' name
)
778 (* get integer parameter
*)
779 val getIParam
= cvt
Int.fromString
780 (* get
real parameter
*)
781 val getRParam
= cvt
Real.fromString
782 (* get
bool parameter
*)
783 val getBParam
= cvt
Bool.fromString
790 * COPYRIGHT (c
) 1993 by AT
&T Bell Laboratories
. See COPYRIGHT file for details
.
792 * Signature for a simple random number generator
.
801 val random
: real -> real
802 (* Given seed
, return value randMin
<= v
<= randMax
803 * Iteratively using the value returned by random
as the
804 * next seed to random will produce a sequence
of pseudo
-random
808 val mkRandom
: real -> unit
-> real
809 (* Given seed
, return function generating a sequence
of
810 * random numbers randMin
<= v
<= randMax
813 val norm
: real -> real
814 (* r
-> r
/ (randMax
+ 1.0) *)
816 val range
: (int * int) -> real -> int
817 (* Map v
, randMin
<= v
<= randMax to integer range
[i
,j
]
825 * COPYRIGHT (c
) 1991 by AT
&T Bell Laboratories
. See COPYRIGHT file for details
827 * Random number generator taken from Paulson
, pp
. 170-171.
828 * Recommended by Stephen K
. Park
and Keith W
. Miller
,
829 * Random number generators
: good ones are hard to find
,
830 * CACM
31 (1988), 1192-1201
832 * Note
: The Random
structure provides a better generator
.
835 structure Rand
: RAND
=
838 (* real number version for systems
with 46-bit mantissas
*)
839 val a
= 16807.0 and m
= 2147483647.0
842 val randMax
= m
- 1.0
844 fun random seed
= let
847 t
- m
* real(floor(t
/m
))
850 fun mkRandom seed
= let
853 fn () => (seed
:= random (!seed
); !seed
)
860 then raise Fail
"Rand.range"
862 val R
= real(j
- i
+ 1)
864 fn r
=> i
+ floor(R
*(r
/m
))
867 val R
= (real j
)-ri
+1.0
869 fn r
=> floor(ri
+ R
*(r
/m
))
876 * COPYRIGHT (c
) 1993, AT
&T Bell Laboratories
.
878 * This is the main driver for the Barnse
-HutN
-body code
.
881 functor Main (V
: VECTOR
) : sig
887 val srand
: int -> unit
888 (* reset the random number generator
*)
890 val testdata
: int -> S
.body list
891 (* generate the Plummer model data
*)
894 output
: {n2bcalc
:int, nbccalc
:int, nstep
:int, selfint
:int, tnow
:real}
896 bodies
: S
.body list
, tnow
: real, tstop
: real,
897 dtime
: real, eps
: real, tol
: real,
898 rmin
: real V
.vec
, rsize
: real
901 val doit
: unit
-> unit
906 structure S
= Space(V
)
907 structure L
= Load(S
)
908 structure G
= Grav(S
)
909 structure DataIO
= DataIO(S
)
911 (* some math utilities
*)
912 (** NOTE
: these are part
of the Math
structure in the new basis
*)
913 val pi
= 3.14159265358979323846
915 if Real.==(y
, 0.0) then 1.0 else Math
.exp (y
* Math
.ln x
)
921 fun srand s
= (seed
:= real s
)
922 fun xrand (xl
, xh
) = let
923 val r
= Rand
.random (! seed
)
926 xl
+ (((xh
- xl
) * r
) / 2147483647.0)
930 (* default parameter values
*)
932 (* file names for input
/output
*)
933 "in=", (* snapshot
of initial conditions
*)
934 "out=", (* stream
of output snapshots
*)
936 (* params
, used
if no input specified
, to make a Plummer Model
*)
937 "nbody=128", (* number
of particles to generate
*)
938 "seed=123", (* random number generator seed
*)
940 (* params to control N
-body integration
*)
941 "dtime=0.025", (* integration time
-step
*)
942 "eps=0.05", (* usual potential softening
*)
943 "tol=1.0", (* cell subdivision tolerence
*)
944 "fcells=0.75", (* cell allocation parameter
*)
946 "tstop=2.0", (* time to stop integration
*)
947 "dtout=0.25", (* data
-output interval
*)
949 "debug=false", (* turn on debugging messages
*)
950 "VERSION=1.0" (* JEB
06 March
1988 *)
953 (* pick a random point on a sphere
of specified radius
. *)
954 fun pickshell rad
= let
956 val vec
= V
.tabulate (fn _
=> xrand(~
1.0, 1.0))
957 val rsq
= V
.dotvp(vec
, vec
)
961 else V
.mulvs (vec
, rad
/ Math
.sqrt(rsq
))
967 (* generate Plummer model initial conditions for test runs
, scaled
968 * to units such that M
= -4E
= G
= 1 (Henon
, Hegge
, etc
).
969 * See Aarseth
, SJ
, Henon
, M
, & Wielen
, R (1974) Astr
& Ap
, 37, 183.
972 val mfrac
= 0.999 (* mass cut off at mfrac
of total
*)
974 val rsc
= (3.0 * pi
) / 16.0
975 val vsc
= Math
.sqrt(1.0 / rsc
)
976 fun mkBodies (0, cmr
, cmv
, l
) = let
977 (* offset bodies by normalized cm coordinates
. Also
, reverse
978 * the list to get the same order
of bodies
as in the C version
.
980 val cmr
= V
.divvs(cmr
, rn
)
981 val cmv
= V
.divvs(cmv
, rn
)
983 |
norm ((p
as S
.Body
{pos
, vel
, ...}) :: r
, l
) = (
984 pos
:= V
.subv(!pos
, cmr
);
985 vel
:= V
.subv(!vel
, cmv
);
990 |
mkBodies (i
, cmr
, cmv
, l
) = let
991 val r
= 1.0 / Math
.sqrt (pow(xrand(0.0, mfrac
), ~
2.0/3.0) - 1.0)
992 val pos
= pickshell (rsc
* r
)
993 fun vN () = let (* von Neumann technique
*)
994 val x
= xrand(0.0,1.0)
995 val y
= xrand(0.0,0.1)
997 if (y
> x
*x
* (pow (1.0-x
*x
, 3.5))) then vN () else x
999 val v
= ((Math
.sqrt
2.0) * vN()) / pow(1.0 + r
*r
, 0.25)
1000 val vel
= pickshell (vsc
* v
)
1009 mkBodies (i
-1, V
.addv(cmr
, pos
), V
.addv(cmv
, vel
), body
:: l
)
1012 mkBodies (n
, V
.zerov
, V
.zerov
, [])
1015 (* startup hierarchical N
-body code
. This either reads
in or generates
1016 * an initial set
of bodies
, and other parameters
.
1019 fun startrun argv
= let
1020 val _
= GetParam
.initParam(argv
, defaults
)
1021 val {nbody
, bodies
, tnow
, headline
} = (case (GetParam
.getParam
"in")
1023 val nbody
= GetParam
.getIParam
"nbody"
1026 then raise Fail
"startrun: absurd nbody"
1028 srand (GetParam
.getIParam
"seed");
1030 bodies
= testdata nbody
,
1032 headline
= "Hack code: Plummer model"
1035 | fname
=> DataIO
.inputData fname
1040 headline
= headline
,
1041 outfile
= GetParam
.getParam
"out",
1042 dtime
= GetParam
.getRParam
"dtime",
1043 eps
= GetParam
.getRParam
"eps",
1044 tol
= GetParam
.getRParam
"tol",
1046 tstop
= GetParam
.getRParam
"tstop",
1047 dtout
= GetParam
.getRParam
"dtout",
1048 debug
= GetParam
.getBParam
"debug",
1049 rmin
= V
.tabulate (fn _
=> ~
2.0),
1054 (* advance N
-body system one time
-step
. *)
1055 fun stepSystem output
{plist
, dtime
, eps
, nstep
, rmin
, rsize
, tnow
, tol
} = let
1056 val dthf
= 0.5 * dtime
1057 val S
.Space
{rmin
, rsize
, root
} = L
.makeTree (plist
, rmin
, rsize
)
1058 (* recalculate accelaration
*)
1059 fun recalc ([], n2bcalc
, nbccalc
, selfint
) = (n2bcalc
, nbccalc
, selfint
)
1060 |
recalc (p
::r
, n2bcalc
, nbccalc
, selfint
) = let
1061 val S
.Body
{acc
as ref acc1
, vel
, ...} = p
1062 val {n2bterm
, nbcterm
, skipSelf
} = G
.hackGrav
{
1063 body
= p
, root
= root
, rsize
= rsize
, eps
= eps
, tol
= tol
1067 then (* use change
in accel to make
2nd order
*)
1068 (* correction to vel
. *)
1069 vel
:= V
.addv(!vel
, V
.mulvs(V
.subv(!acc
, acc1
), dthf
))
1071 recalc (r
, n2bcalc
+n2bterm
, nbccalc
+nbcterm
,
1072 if skipSelf
then selfint
else (selfint
+1))
1074 (* advance bodies
*)
1075 fun advance (S
.Body
{pos
, acc
, vel
, ...}) = let
1076 val dvel
= V
.mulvs (!acc
, dthf
)
1077 val vel1
= V
.addv (!vel
, dvel
)
1078 val dpos
= V
.mulvs (vel1
, dtime
)
1080 pos
:= V
.addv (!pos
, dpos
);
1081 vel
:= V
.addv (vel1
, dvel
)
1083 val (n2bcalc
, nbccalc
, selfint
) = recalc (plist
, 0, 0, 0)
1085 output
{nstep
=nstep
, tnow
=tnow
, n2bcalc
=n2bcalc
, nbccalc
=nbccalc
, selfint
=selfint
};
1087 (nstep
+1, tnow
+ dtime
)
1090 (* given an initial configuration
, run the simulation
*)
1092 output
, bodies
, tnow
, tstop
,
1093 dtime
, eps
, tol
, rsize
, rmin
1095 val step
= stepSystem output
1096 fun loop (nstep
, tnow
) = if (tnow
< tstop
+ (0.1 * dtime
))
1098 plist
= bodies
, dtime
= dtime
, eps
= eps
, nstep
= nstep
,
1099 rmin
= rmin
, rsize
= rsize
, tnow
= tnow
, tol
= tol
1107 val { nbody
, bodies
, headline
, outfile
,
1108 dtime
, eps
, tol
, tnow
, tstop
, dtout
,
1111 fun output
{nstep
, tnow
, n2bcalc
, nbccalc
, selfint
} = DataIO
.output
{
1112 bodies
= bodies
, nbody
= nbody
,
1113 n2bcalc
= n2bcalc
, nbccalc
= nbccalc
,
1114 selfint
= selfint
, tnow
= tnow
1118 outfile
= outfile
, headline
= headline
, nbody
= nbody
, tnow
= tnow
,
1119 dtime
= dtime
, eps
= eps
, tol
= tol
, dtout
= dtout
, tstop
= tstop
1122 output
=output
, bodies
=bodies
, tnow
=tnow
, tstop
=tstop
,
1123 dtime
=dtime
, eps
=eps
, tol
=tol
, rsize
=rsize
, rmin
=rmin
1131 * COPYRIGHT (c
) 1993, AT
&T Bell Laboratories
.
1133 * 3 dimensional vector arithmetic
.
1136 structure Vector3
: VECTOR
=
1139 type 'a vec
= {x
: 'a
, y
: 'a
, z
: 'a
}
1140 type realvec
= real vec
1144 fun tabulate f
= {x
= f
0, y
= f
1, z
= f
2}
1146 val zerov
= {x
= 0.0, y
= 0.0, z
= 0.0}
1147 fun equal({x
, y
, z
}, {x
=x1
, y
=y1
, z
=z1
}) =
1148 Real.==(x
, x1
) andalso Real.==(y
, y1
) andalso Real.==(z
, z1
)
1150 fun addv ({x
=x1
, y
=y1
, z
=z1
} : realvec
, {x
=x2
, y
=y2
, z
=z2
}) =
1151 {x
=x1
+x2
, y
=y1
+y2
, z
=z1
+z2
}
1153 fun subv ({x
=x1
, y
=y1
, z
=z1
} : realvec
, {x
=x2
, y
=y2
, z
=z2
}) =
1154 {x
=x1
-x2
, y
=y1
-y2
, z
=z1
-z2
}
1156 fun dotvp ({x
=x1
, y
=y1
, z
=z1
} : realvec
, {x
=x2
, y
=y2
, z
=z2
}) =
1157 x1
*x2
+ y1
*y2
+ z1
*z2
1159 fun crossvp ({x
=x1
, y
=y1
, z
=z1
} : realvec
, {x
=x2
, y
=y2
, z
=z2
}) =
1160 {x
= y1
*z2
- z1
*y2
, y
= x1
*z2
- z1
*x2
, z
= x1
*y2
- y1
*x2
}
1162 fun addvs ({x
, y
, z
} : realvec
, s
) = {x
=x
+s
, y
=y
+s
, z
=z
+s
}
1164 fun mulvs ({x
, y
, z
} : realvec
, s
) = {x
=x
*s
, y
=y
*s
, z
=z
*s
}
1166 fun divvs ({x
, y
, z
} : realvec
, s
) = {x
=x
/s
, y
=y
/s
, z
=z
/s
}
1168 fun mapv f
{x
, y
, z
} = {x
= f x
, y
= f y
, z
= f z
}
1170 fun map3v
f ({x
=x1
, y
=y1
, z
=z1
}, {x
=x2
, y
=y2
, z
=z2
}, {x
=x3
, y
=y3
, z
=z3
}) =
1171 {x
= f(x1
, x2
, x3
), y
= f(y1
, y2
, y3
), z
= f(z1
, z2
, z3
)}
1173 fun foldv f
{x
, y
, z
} init
= f(z
, f(y
, f(x
, init
)))
1175 fun format
{lp
, rp
, sep
, cvt
} {x
, y
, z
} = String.concat
[
1176 lp
, cvt x
, sep
, cvt y
, sep
, cvt z
, rp
1179 fun explode
{x
, y
, z
} = [x
, y
, z
]
1181 fun implode
[x
, y
, z
] = {x
=x
, y
=y
, z
=z
}
1182 | implode _
= raise Fail
"implode: bad dimension"
1185 m00
: real, m01
: real, m02
: real,
1186 m10
: real, m11
: real, m12
: real,
1187 m20
: real, m21
: real, m22
: real
1191 m00
= 0.0, m01
= 0.0, m02
= 0.0,
1192 m10
= 0.0, m11
= 0.0, m12
= 0.0,
1193 m20
= 0.0, m21
= 0.0, m22
= 0.0
1196 fun addm (a
: matrix
, b
: matrix
) = {
1197 m00
=(#m00 a
+ #m00 b
), m01
=(#m01 a
+ #m01 b
), m02
=(#m02 a
+ #m02 b
),
1198 m10
=(#m10 a
+ #m10 b
), m11
=(#m11 a
+ #m11 b
), m12
=(#m12 a
+ #m12 b
),
1199 m20
=(#m20 a
+ #m20 b
), m21
=(#m21 a
+ #m21 b
), m22
=(#m22 a
+ #m22 b
)
1202 fun outvp ({x
=a0
, y
=a1
, z
=a2
} : realvec
, {x
=b0
, y
=b1
, z
=b2
}) = {
1203 m00
=(a0
*b0
), m01
=(a0
*b1
), m02
=(a0
*b2
),
1204 m10
=(a1
*b0
), m11
=(a1
*b1
), m12
=(a1
*b2
),
1205 m20
=(a2
*b0
), m21
=(a2
*b1
), m22
=(a2
*b2
)
1212 val doit
: int -> unit
1213 val testit
: TextIO.outstream
-> unit
1215 (* load file for bmark version
*)
1231 structure Main
: BMARK
=
1233 structure M3
= Main(Vector3
);
1235 val name
= "Barnes-Hut (3d)"
1237 fun testit strm
= ()
1242 output
= fn _
=> (),
1243 bodies
= M3
.testdata n
,
1244 tnow
= 0.0, tstop
= 2.0,
1245 dtime
= 0.025, eps
= 0.05, tol
= 1.0,
1246 rmin
= M3
.S
.V
.tabulate (fn _
=> ~
2.0),