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

Annotation of /sml/trunk/benchmarks/todo/bh/bh.sml

Parent Directory Parent Directory | Revision Log Revision Log


Revision 193 - (view) (download)

1 : monnier 193 (* hierachical N-body code, adapted from Barnes-Hut C version *)
2 :    
3 :     structure BarnesHut =
4 :     struct
5 :     open bhMath bhUtil Array
6 :    
7 :     val IMAX = 536870912 (* 2^29 *)
8 :     val IMAXrs1 = Bits.rshift(IMAX,1)
9 :     val NDIM = 3
10 :     val NSUB = Bits.lshift(1,NDIM)
11 :    
12 :     datatype content = Body of {vel:real vec ref,
13 :     acc:real vec ref,
14 :     phi:real ref} |
15 :     Cell of node array
16 :     and
17 :     node = Empty |
18 :     Node of {mass:real,
19 :     pos: real vec ref,
20 :     contents:content}
21 :    
22 :     fun putCell (c as Cell x) k v = (update(x,k,v);
23 :     c)
24 :     fun getCell (Cell x) k = sub(x,k)
25 :    
26 :     fun pickshell rad =
27 :     let fun pickvec () =
28 :     let val rv = randvec (~1.0,1.0)
29 :     in
30 :     if (dotvp rv rv) > 1.0 then
31 :     pickvec ()
32 :     else
33 :     rv
34 :     end
35 :     val vec = pickvec ()
36 :     in
37 :     mulvs vec (rad / (sqrt (dotvp vec vec)))
38 :     end
39 :    
40 :     fun buildtest n = (* builds Plummer model *)
41 :     let val rn = real n
42 :     val bodymass = 1.0 / (rn)
43 :     val mfrac = 0.9999
44 :     val sqrt2 = sqrt 2.0
45 :     val rsc = 3.0 * PI / 16.0
46 :     val vsc = sqrt (1.0 / rsc)
47 :     fun vN () = (* von Neumann technique *)
48 :     let val x = xrand(0.0,1.0)
49 :     val y = xrand(0.0,0.1)
50 :     in
51 :     if (y > x*x * (pow (1.0-x*x) 3.5)) then vN ()
52 :     else x
53 :     end
54 :     fun mkbody _ =
55 :     let val r = 1.0/sqrt(pow(xrand(0.00001,mfrac))(~2.0/3.0)-1.0)
56 :     in
57 :     Node{mass=bodymass,
58 :     pos=ref (pickshell (rsc * r)),
59 :     contents=Body{vel=ref
60 :     (pickshell (vsc * (sqrt2 * (vN ()) /
61 :     (pow (1.0+r*r) 0.25)))),
62 :     acc=ref realzerovec,
63 :     phi=ref 0.0}}
64 :     end
65 :     fun scale (n as Node{pos,contents=Body{vel,...},...}) ps vs =
66 :     (pos := subv (!pos) ps;
67 :     vel := subv (!vel) vs;
68 :     n)
69 :     val bodies = map mkbody (intsto n)
70 :     val cmr = fold (fn (Node{pos,...},v) => addv (!pos) v)
71 :     bodies realzerovec
72 :     val cmr' = divvs (divvs cmr rn) rn
73 :     val cmv = fold (fn (Node{contents=Body{vel,...},...},v) =>
74 :     addv (!vel) v)
75 :     bodies realzerovec
76 :     val cmv' = divvs (divvs cmv rn) rn
77 :     in
78 :     map (fn b => scale b cmr' cmv') bodies
79 :     end
80 :    
81 :     fun intcoord rp rmin rsize =
82 :     let val xsc = divvs (subv rp rmin) rsize
83 :     exception notintcoord
84 :     in
85 :     SOME (map (fn x => if (x >= 0.0 andalso x < 1.0) then
86 :     floor (real IMAX * x)
87 :     else
88 :     raise notintcoord)
89 :     xsc)
90 :     handle notintcoord => NONE
91 :     | x => raise x
92 :     end
93 :    
94 :     fun cellindex iv l =
95 :     let fun aux [] k = 0
96 :     | aux (v::vs) k = (if (Bits.andb(l,v)) <> 0 then
97 :     Bits.rshift(NSUB,k+1)
98 :     else
99 :     0) + (aux vs (k+1))
100 :     in
101 :     aux iv 0
102 :     end
103 :    
104 :     fun expandbox (nd as Node{pos,...}) (state as (rmin,rsize,root)) =
105 :     let val ic = intcoord (!pos) rmin rsize
106 :     in
107 :     case ic of
108 :     NONE =>
109 :     let val rmid = addvs rmin (0.5 * rsize)
110 :     val rmin' = mapvec3 (fn (x,y,z) => if x < y then
111 :     z - rsize
112 :     else
113 :     z)
114 :     (!pos) rmid rmin
115 :     val rsize' = 2.0 * rsize
116 :     fun mksub v r =
117 :     let val SOME x = (intcoord v rmin' rsize')
118 :     val k = cellindex x IMAXrs1
119 :     val a = Cell (array(NSUB,Empty))
120 :     in
121 :     putCell a k r
122 :     end
123 :     in
124 :     expandbox nd (case root of
125 :     Empty =>
126 :     (rmin',rsize',root)
127 :     | Node{...} =>
128 :     (rmin',
129 :     rsize',
130 :     Node{mass=0.0,
131 :     pos=ref realzerovec,
132 :     contents=mksub rmid root}))
133 :     end
134 :     | SOME _ => state
135 :     end
136 :    
137 :     fun walk tree (acc0,eps,phi0,pos0,pskip,rsize,tol) =
138 :     let fun subdivp (Node{contents=Body _,...}) _ _ = false
139 :     | subdivp (Node{pos=pos,...}) dsq (pos0,tol) =
140 :     let val tolsq = tol * tol
141 :     val dr = subv (!pos) pos0
142 :     val drsq = dotvp dr dr
143 :     in
144 :     (tolsq * drsq) < dsq
145 :     end
146 :     fun gravsub (Node{pos=pos,mass=mass,...}) (pos0,eps,acc0,phi0) =
147 :     let val dr = subv (!pos) pos0
148 :     val drsq = (dotvp dr dr) + eps * eps
149 :     val drabs = sqrt drsq
150 :     val phii = mass / drabs
151 :     val phi0' = phi0 - phii
152 :     val mor3 = phii / drsq
153 :     val acc0' = addv acc0 (mulvs dr mor3)
154 :     in
155 :     (acc0',phi0')
156 :     end
157 :     fun walksub p dsq (state as (acc0,eps,phi0,pos0,pskip,tol)) =
158 :     if subdivp p dsq (pos0,tol) then
159 :     let val Node{contents=Cell subcells,...} = p
160 :     in
161 :     foldarray (fn (n as Node _,s) =>
162 :     walksub n (dsq / 4.0) s
163 :     | (_,s) => s)
164 :     subcells
165 :     state
166 :     end
167 :     else if p = pskip then state
168 :     else
169 :     let val (acc0',phi0') = gravsub p (pos0,eps,acc0,phi0)
170 :     in
171 :     (acc0',eps,phi0',pos0,pskip,tol)
172 :     end
173 :     val (acc0',_,phi0',_,_,_) = walksub
174 :     tree (rsize * rsize)
175 :     (acc0,eps,phi0,pos0,pskip,tol)
176 :     in
177 :     (acc0',phi0')
178 :     end
179 :    
180 :     fun grav (p as Node{pos=pos,contents=Body{phi=phi,acc=acc,...},...})
181 :     (root,eps,rsize,tol) =
182 :     let val (acc0',phi0') = walk root
183 :     (realzerovec,eps,0.0,!pos,p,rsize,tol)
184 :     in
185 :     phi := phi0';
186 :     acc := acc0'
187 :     end
188 :    
189 :     fun loadtree (n as Node{pos=posn,...}) (rmin,rsize,root) =
190 :     let val SOME xp = (intcoord (!posn) rmin rsize)
191 :     val l = IMAXrs1
192 :     fun aux Empty _ = n
193 :     | aux (n as Node{contents=Body{...},pos=posb,...}) l =
194 :     let val SOME xq = (intcoord (!posb) rmin rsize)
195 :     val k = cellindex xq l
196 :     val a = Cell (array(NSUB,Empty))
197 :     in
198 :     aux (Node{mass=0.0,
199 :     pos=ref realzerovec,
200 :     contents=putCell a k n})
201 :     l
202 :     end
203 :     | aux (n as Node{contents=c as Cell x,...}) l =
204 :     let val i = cellindex xp l
205 :     in
206 :     putCell c i (aux (getCell c i)
207 :     (Bits.rshift(l,1)));
208 :     n
209 :     end
210 :     in
211 :     (rmin,rsize,aux root l)
212 :     end
213 :    
214 :     fun cellmass (Node{contents=Cell subcells,...}) =
215 :     let val () = maparray cellmass subcells
216 :     val m = foldarray (fn (Node{mass=mass,...},x) => mass+x
217 :     | (_,x) => x)
218 :     subcells
219 :     0.0
220 :     in
221 :     Node{mass=m,pos=ref realzerovec,contents=Cell subcells}
222 :     end
223 :     | cellmass x = x
224 :    
225 :     fun cellmoment (nd as Node{contents=Cell cells,pos=pos,mass=mass,...}) =
226 :     let val () = maparray cellmoment cells
227 :     val moment = foldarray (fn (Node{pos=p,mass=m,...},v) =>
228 :     addv v (mulvs (!p) m)
229 :     | (_,v) => v)
230 :     cells
231 :     (!pos)
232 :     in
233 :     pos := divvs moment mass;
234 :     nd
235 :     end
236 :     | cellmoment x = x
237 :    
238 :    
239 :     fun maketree nodelist x =
240 :     let fun aux (n as Node{mass,...},state) =
241 :     if mass = 0.0 then
242 :     state
243 :     else
244 :     loadtree n (expandbox n state)
245 :     in
246 :     fold aux nodelist x
247 :     end
248 :    
249 :     fun stepsystem plist (dtime,eps,nstep,rmin,rsize,tnow,tol) =
250 :     let val (rmin,rsize,tree) = maketree plist (rmin,rsize,Empty)
251 :     val tree' = cellmoment (cellmass tree)
252 :     val dthf = 0.5 * dtime
253 :     fun recalc (n as Node{contents=Body{acc=acc,vel=vel,...},...}) =
254 :     let val acc1 = !acc
255 :     val _ = grav n (tree',eps,rsize,tol)
256 :     in
257 :     if nstep > 0 then
258 :     let val dacc = subv (!acc) acc1
259 :     val dvel = mulvs dacc dthf
260 :     in
261 :     vel := addv (!vel) dvel
262 :     end
263 :     else
264 :     ()
265 :     end
266 :     fun adv (Node{pos=pos,contents=Body{acc=acc,vel=vel,...},...}) =
267 :     let val dvel = mulvs (!acc) dthf
268 :     val vel1 = addv (!vel) dvel
269 :     val dpos = mulvs vel1 dtime
270 :     in
271 :     pos := addv (!pos) dpos;
272 :     vel := addv vel1 dvel
273 :     end
274 :     in
275 :     app recalc plist;
276 :     app adv plist;
277 :     (nstep+1,tnow+dtime)
278 :     end
279 :    
280 :     fun dump (Node{mass=m,pos=p,contents=Body{vel=v,...},...}) =
281 :     (print "Mass=";
282 :     print m;
283 :     print "\n";
284 :     print "pos=";
285 :     realprvec (!p);
286 :     print "vel=";
287 :     realprvec (!v))
288 :    
289 :     fun testit n = (* test a model for n bodies *)
290 :     let val bodies = buildtest n
291 :     val dtime = 0.025
292 :     val eps = 0.05
293 :     val nstep = 0
294 :     val tol = 1.0
295 :     val tstop = (* 0.5 *) 2.0
296 :     val tnow = 0.0
297 :     val rmin = mkvs ~2.0
298 :     val rsize = ~2.0 * getx rmin
299 :     fun doit (state as (dtime,eps,nstep,rmin,rsize,tnow,tol)) =
300 :     if tnow < tstop + 0.1 * dtime then
301 :     let fun status step time =
302 :     (print "stepping system at time ";
303 :     print time;
304 :     print " to next step ";
305 :     print step;
306 :     print "\n")
307 :     val _ = status nstep tnow
308 :     (* val _ = app dump bodies *)
309 :     val (nstep',tnow') = stepsystem bodies state
310 :     in
311 :     doit (dtime,eps,nstep',rmin,rsize,tnow',tol)
312 :     end
313 :     else
314 :     ()
315 :     in
316 :     doit (dtime,eps,nstep,rmin,rsize,tnow,tol)
317 :     end
318 :     end
319 :    

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