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

Annotation of /sml/trunk/benchmarks/programs/simple/simple.sml

Parent Directory Parent Directory | Revision Log Revision Log


Revision 193 - (view) (download)

1 : monnier 193 (* Simple
2 :     * error: grid_max < 5
3 :     *)
4 :     functor Simple(val grid_max: int val step_count: int) : BMARK =
5 :     struct
6 :    
7 :     fun fold f [] = (fn b => b)
8 :     | fold f (a::r) = (fn b => let fun f2(e,[]) = f(e,b)
9 :     | f2(e,a::r) = f(e,f2(a,r))
10 :     in f2(a,r)
11 :     end)
12 :    
13 :    
14 :     fun min (x:real,y:real) = if x<y then x else y
15 :     fun max (x:real,y:real) = if x<y then y else x
16 :     exception MaxList
17 :     exception MinList
18 :     exception SumList
19 :     fun max_list [] = raise MaxList | max_list l = fold max l (hd l)
20 :     fun min_list [] = raise MinList | min_list l = fold min l (hd l)
21 :     fun sum_list [] = raise SumList
22 :     | sum_list (l:real list) = fold (op +) l 0.0
23 :    
24 :     fun for {from=start:int,step=delta:int, to=endd:int} body =
25 :     if delta>0 andalso endd>=start then
26 :     let fun f x = if x > endd then () else (body x; f(x+delta))
27 :     in f start
28 :     end
29 :     else if endd<=start then
30 :     let fun f x = if x < endd then () else (body x; f(x+delta))
31 :     in f start
32 :     end
33 :     else ()
34 :     fun from(n,m) = if n>m then [] else n::from(n+1,m)
35 :     fun flatten [] = []
36 :     | flatten (x::xs) = x @ flatten xs
37 :     fun pow(x:real,y:int) = if y = 0 then 1.0 else x * pow(x,y-1)
38 :     fun array2(bounds as ((l1,u1),(l2,u2)),v) =
39 :     (Array2.array((u1-l1+1, u2-l2+1),v), bounds)
40 :     fun sub2((A,((lb1:int,ub1:int),(lb2:int,ub2:int))),(k,l)) =
41 :     Array2.sub(A, (k-lb1, l-lb2))
42 :     fun update2((A,((lb1,_),(lb2,_))),(k,l), v) = Array2.update(A,(k-lb1,l-lb2),v)
43 :     fun bounds2(_,b) = b
44 :     fun printarray2 (A as (M:real Array2.array2,((l1,u1),(l2,u2)))) =
45 :     for {from=l1,step=1,to=u1} (fn i =>
46 :     (print "[";
47 :     for {from=l2,step=1,to=u2-1} (fn j =>
48 :     print (Real.toString (sub2(A,(i,j))) ^ ", "));
49 :     print (Real.toString (sub2(A,(i,u2))) ^ "]\n")))
50 :     fun array1((l,u),v) = (Array.array(u-l+1,v),(l,u))
51 :     fun sub1((A,(l:int,u:int)),i:int) = Array.sub(A,i-l)
52 :     fun update1((A,(l,_)),i,v) = Array.update(A,i-l,v)
53 :     fun bounds1(_,b) = b
54 :    
55 :     (*
56 :     * Specification of the state variable computation
57 :     *)
58 :     val grid_size = ((2,grid_max), (2,grid_max))
59 :    
60 :     fun north (k,l) = (k-1,l)
61 :     fun south (k,l) = (k+1,l)
62 :    
63 :     fun east (k,l) = (k,l+1)
64 :     fun west (k,l) = (k,l-1)
65 :    
66 :     val northeast = north o east
67 :     val southeast = south o east
68 :     val northwest = north o west
69 :     val southwest = south o west
70 :    
71 :     fun farnorth x = (north o north ) x
72 :     fun farsouth x = (south o south) x
73 :     fun fareast x = (east o east) x
74 :     fun farwest x = (west o west) x
75 :    
76 :     fun zone_A(k,l) = (k,l)
77 :     fun zone_B(k,l) = (k+1,l)
78 :    
79 :     fun zone_C(k,l) = (k+1,l+1)
80 :     fun zone_D(k,l) = (k,l+1)
81 :    
82 :     val zone_corner_northeast = north
83 :     val zone_corner_northwest = northwest
84 :     fun zone_corner_southeast zone = zone
85 :     val zone_corner_southwest = west
86 :    
87 :     val ((kmin,kmax),(lmin,lmax)) = grid_size
88 :     val dimension_all_nodes = ((kmin-1,kmax+1),(lmin-1,lmax+1))
89 :     fun for_all_nodes f =
90 :     for {from=kmin-1, step=1, to=kmax+1} (fn k =>
91 :     for {from=lmin-1, step=1, to=lmax+1} (fn l => f k l))
92 :    
93 :     val dimension_interior_nodes = ((kmin,kmax),(lmin,lmax))
94 :     fun for_interior_nodes f =
95 :     for {from=kmin, step=1, to=kmax} (fn k =>
96 :     for {from=lmin, step=1, to=lmax} (fn l => f k l))
97 :    
98 :     val dimension_all_zones = ((kmin,kmax+1),(lmin,lmax+1))
99 :     fun for_all_zones f =
100 :     for {from=kmin, step=1, to=kmax+1} (fn k =>
101 :     for {from=lmin, step=1, to=lmax+1} (fn l => f (k,l)))
102 :    
103 :     val dimension_interior_zones = ((kmin+1,kmax),(lmin+1,lmax))
104 :     fun for_interior_zones f =
105 :     for {from=kmin+1, step=1, to=kmax} (fn k =>
106 :     for {from=lmin+1, step=1, to=lmax} (fn l => f (k,l)))
107 :    
108 :     fun map_interior_nodes f =
109 :     flatten(map (fn k => (map (fn l => f (k,l))
110 :     (from(lmin,lmax))))
111 :     (from(kmin,kmax)))
112 :     fun map_interior_zones f =
113 :     flatten(map (fn k => (map (fn l => f (k,l))
114 :     (from(lmin+1,lmax))))
115 :     (from(kmin+1,kmax)))
116 :    
117 :     fun for_north_ward_interior_zones f =
118 :     for {from=kmax, step= ~1, to=kmin+1} (fn k =>
119 :     for {from=lmin+1, step=1, to=lmax} (fn l => f (k,l)))
120 :     fun for_west_ward_interior_zones f =
121 :     for {from=kmin+1, step=1, to=kmax} (fn k =>
122 :     for {from=lmax, step= ~1, to=lmin+1} (fn l => f (k,l)))
123 :    
124 :    
125 :     fun for_north_zones f = for {from=lmin, step=1, to=lmax+1} (fn l => f (kmin,l))
126 :     fun for_south_zones f = for {from=lmin+1, step=1, to=lmax} (fn l => f (kmax+1,l))
127 :     fun for_east_zones f = for {from=kmin+1, step=1, to=kmax+1}(fn k => f (k,lmax+1))
128 :     fun for_west_zones f = for {from=kmin+1, step=1, to=kmax+1}(fn k => f (k,lmin))
129 :    
130 :     fun reflect dir node A = sub2(A, dir node)
131 :     val reflect_north = fn x => reflect north x
132 :     val reflect_south = fn x => reflect south x
133 :     val reflect_east = fn x => reflect east x
134 :     val reflect_west = fn x => reflect west x
135 :    
136 :     fun for_north_nodes f =
137 :     for {from=lmin, step=1, to=lmax-1} (fn l => f (kmin-1,l))
138 :     fun for_south_nodes f =
139 :     for {from=lmin, step=1, to=lmax-1} (fn l => f (kmax+1,l))
140 :     fun for_east_nodes f =
141 :     for {from=kmin, step=1, to=kmax-1} (fn k => f (k,lmax+1))
142 :     fun for_west_nodes f =
143 :     for {from=kmin, step=1, to=kmax-1} (fn k => f (k,lmin-1))
144 :    
145 :     val north_east_corner = (kmin-1,lmax+1)
146 :     val north_west_corner = (kmin-1,lmin-1)
147 :     val south_east_corner = (kmax+1,lmax+1)
148 :     val south_west_corner = (kmax+1,lmin-1)
149 :    
150 :     val west_of_north_east = (kmin-1, lmax)
151 :     val west_of_south_east = (kmax+1, lmax)
152 :     val north_of_south_east = (kmax, lmax+1)
153 :     val north_of_south_west = (kmax, lmin-1)
154 :    
155 :    
156 :    
157 :     (*
158 :     * Initialization of parameters
159 :     *)
160 :     val constant_heat_source = 0.0
161 :     val deltat_maximum = 0.01
162 :     val specific_heat = 0.1
163 :     val p_coeffs = let val M = array2(((0,2),(0,2)), 0.0)
164 :     in update2(M, (1,1), 0.06698); M
165 :     end
166 :     val e_coeffs = let val M = array2(((0,2),(0,2)), 0.0)
167 :     in update2(M, (0,1), 0.1); M
168 :     end
169 :     val p_poly = array2(((1,4),(1,5)),p_coeffs)
170 :    
171 :     val e_poly = array2(((1,4),(1,5)), e_coeffs)
172 :    
173 :     val rho_table = let val V = array1((1,3), 0.0)
174 :     in update1(V,2,1.0);
175 :     update1(V,3,100.0);
176 :     V
177 :     end
178 :     val theta_table = let val V = array1((1,4), 0.0)
179 :     in update1(V,2,3.0);
180 :     update1(V,3,300.0);
181 :     update1(V,4,3000.0);
182 :     V
183 :     end
184 :    
185 :     val extract_energy_tables_from_constants = (e_poly,2,rho_table,theta_table)
186 :     val extract_pressure_tables_from_constants = (p_poly,2,rho_table,theta_table)
187 :    
188 :     val nbc = let val M = array2(dimension_all_zones, 1)
189 :     in for {from=lmin+1,step=1,to=lmax} (fn j => update2(M,(kmax+1, j),2));
190 :     update2(M,(kmin,lmin),4);
191 :     update2(M,(kmin,lmax+1),4);
192 :     update2(M,(kmax+1,lmin),4);
193 :     update2(M,(kmax+1,lmax+1),4);
194 :     M
195 :     end
196 :     val pbb = let val A = array1((1,4), 0.0)
197 :     in update1(A,2,6.0); A
198 :     end
199 :     val pb = let val A = array1((1,4), 1.0)
200 :     in update1(A,2,0.0); update1(A,3,0.0); A
201 :     end
202 :     val qb = pb
203 :    
204 :     val all_zero_nodes = array2(dimension_all_nodes, 0.0)
205 :    
206 :     val all_zero_zones = array2(dimension_all_zones, 0.0)
207 :    
208 :    
209 :     (*
210 :     * Positional Coordinates. (page 9-10)
211 :     *)
212 :    
213 :     fun make_position_matrix interior_function =
214 :     let val r' = array2(dimension_all_nodes, 0.0)
215 :     val z' = array2(dimension_all_nodes, 0.0)
216 :     fun boundary_position (rx,zx,ry,zy,ra,za) =
217 :     let val (rax, zax) = (ra - rx, za - zx)
218 :     val (ryx, zyx) = (ry - rx, zy - zx)
219 :     val omega = 2.0*(rax*ryx + zax*zyx)/(ryx*ryx + zyx*zyx)
220 :     val rb = rx - rax + omega*ryx
221 :     val zb = zx - zax + omega*zyx
222 :     in (rb, zb)
223 :     end
224 :    
225 :     fun reflect_node (x_dir, y_dir, a_dir, node) =
226 :     let val rx = reflect x_dir node r'
227 :     val zx = reflect x_dir node z'
228 :     val ry = reflect y_dir node r'
229 :     val zy = reflect y_dir node z'
230 :     val ra = reflect a_dir node r'
231 :     val za = reflect a_dir node z'
232 :     in boundary_position (rx, zx, ry, zy, ra, za)
233 :     end
234 :     fun u2 (rv,zv) n = (update2(r',n,rv); update2(z',n,zv))
235 :     in
236 :     for_interior_nodes (fn k => fn l => u2 (interior_function (k,l)) (k,l));
237 :     for_north_nodes(fn n => u2 (reflect_node(south,southeast,farsouth,n)) n);
238 :     for_south_nodes (fn n => u2(reflect_node(north,northeast,farnorth,n)) n);
239 :     for_east_nodes (fn n => u2(reflect_node(west, southwest, farwest, n)) n);
240 :     for_west_nodes (fn n => u2(reflect_node(east, southeast, fareast, n)) n);
241 :     u2 (reflect_node(south, southwest, farsouth, west_of_north_east))
242 :     west_of_north_east;
243 :     u2 (reflect_node(north, northwest, farnorth, west_of_south_east))
244 :     west_of_south_east;
245 :     u2 (reflect_node(west, northwest, farwest, north_of_south_east))
246 :     north_of_south_east;
247 :     u2 (reflect_node(east, northeast, fareast, north_of_south_west))
248 :     north_of_south_west;
249 :     u2 (reflect_node(southwest, west, farwest, north_east_corner))
250 :     north_east_corner;
251 :     u2 (reflect_node(northwest, west, farwest, south_east_corner))
252 :     south_east_corner;
253 :     u2 (reflect_node(southeast, south, farsouth, north_west_corner))
254 :     north_west_corner;
255 :     u2 (reflect_node(northeast, east, fareast, south_west_corner))
256 :     south_west_corner;
257 :     (r',z')
258 :     end
259 :    
260 :    
261 :    
262 :     (*
263 :     * Physical Properties of a Zone (page 10)
264 :     *)
265 :     fun zone_area_vol ((r,z), zone) =
266 :     let val (r1,z1)=(sub2(r,zone_corner_southwest zone),
267 :     sub2(z,zone_corner_southwest zone))
268 :     val (r2,z2)=(sub2(r,zone_corner_southeast zone),
269 :     sub2(z,zone_corner_southeast zone))
270 :     val (r3,z3)=(sub2(r,zone_corner_northeast zone),
271 :     sub2(z,zone_corner_northeast zone))
272 :     val (r4,z4)=(sub2(r,zone_corner_northwest zone),
273 :     sub2(z,zone_corner_northwest zone))
274 :     val area1 = (r2-r1)*(z3-z1) - (r3-r2)*(z3-z2)
275 :     val radius1 = 0.3333 *(r1+r2+r3)
276 :     val volume1 = area1 * radius1
277 :     val area2 = (r3-r1)*(z4-z3) - (r4-r3)*(z3-z1)
278 :     val radius2 = 0.3333 *(r1+r3+r4)
279 :     val volume2 = area2 * radius2
280 :     in (area1+area2, volume1+volume2)
281 :     end
282 :    
283 :     (*
284 :     * Velocity (page 8)
285 :     *)
286 :     fun make_velocity((u,w),(r,z),p,q,alpha,rho,delta_t) =
287 :     let fun line_integral (p,z,node) : real =
288 :     sub2(p,zone_A node)*(sub2(z,west node) - sub2(z,north node)) +
289 :     sub2(p,zone_B node)*(sub2(z,south node) - sub2(z,west node)) +
290 :     sub2(p,zone_C node)*(sub2(z,east node) - sub2(z,south node)) +
291 :     sub2(p,zone_D node)*(sub2(z,north node) - sub2(z,east node))
292 :     fun regional_mass node =
293 :     0.5 * (sub2(rho, zone_A node)*sub2(alpha,zone_A node) +
294 :     sub2(rho, zone_B node)*sub2(alpha,zone_B node) +
295 :     sub2(rho, zone_C node)*sub2(alpha,zone_C node) +
296 :     sub2(rho, zone_D node)*sub2(alpha,zone_D node))
297 :     fun velocity node =
298 :     let val d = regional_mass node
299 :     val n1 = ~(line_integral(p,z,node)) - line_integral(q,z,node)
300 :     val n2 = line_integral(p,r,node) + line_integral(q,r,node)
301 :     val u_dot = n1/d
302 :     val w_dot = n2/d
303 :     in (sub2(u,node)+delta_t*u_dot, sub2(w,node)+delta_t*w_dot)
304 :     end
305 :     val U = array2(dimension_interior_nodes,0.0)
306 :     val W = array2(dimension_interior_nodes,0.0)
307 :     in for_interior_nodes (fn k => fn l => let val (uv,wv) = velocity (k,l)
308 :     in update2(U,(k,l),uv);
309 :     update2(W,(k,l),wv)
310 :     end);
311 :     (U,W)
312 :     end
313 :    
314 :    
315 :    
316 :     fun make_position ((r,z),delta_t,(u',w')) =
317 :     let fun interior_position node =
318 :     (sub2(r,node) + delta_t*sub2(u',node),
319 :     sub2(z,node) + delta_t*sub2(w',node))
320 :     in make_position_matrix interior_position
321 :     end
322 :    
323 :    
324 :     fun make_area_density_volume(rho, s, x') =
325 :     let val alpha' = array2(dimension_all_zones, 0.0)
326 :     val s' = array2(dimension_all_zones, 0.0)
327 :     val rho' = array2(dimension_all_zones, 0.0)
328 :     fun interior_area zone =
329 :     let val (area, vol) = zone_area_vol (x', zone)
330 :     val density = sub2(rho,zone)*sub2(s,zone) / vol
331 :     in (area,vol,density)
332 :     end
333 :     fun reflect_area_vol_density reflect_function =
334 :     (reflect_function alpha',reflect_function s',reflect_function rho')
335 :     fun update_asr (zone,(a,s,r)) = (update2(alpha',zone,a);
336 :     update2(s',zone,s);
337 :     update2(rho',zone,r))
338 :     fun r_area_vol_den (reflect_dir,zone) =
339 :     let val asr = reflect_area_vol_density (reflect_dir zone)
340 :     in update_asr(zone, asr)
341 :     end
342 :     in
343 :     for_interior_zones (fn zone => update_asr(zone, interior_area zone));
344 :     for_south_zones (fn zone => r_area_vol_den(reflect_north, zone));
345 :     for_east_zones (fn zone => r_area_vol_den(reflect_west, zone));
346 :     for_west_zones (fn zone => r_area_vol_den(reflect_east, zone));
347 :     for_north_zones (fn zone => r_area_vol_den(reflect_south, zone));
348 :     (alpha', rho', s')
349 :     end
350 :    
351 :    
352 :     (*
353 :     * Artifical Viscosity (page 11)
354 :     *)
355 :     fun make_viscosity(p,(u',w'),(r',z'), alpha',rho') =
356 :     let fun interior_viscosity zone =
357 :     let fun upper_del f =
358 :     0.5 * ((sub2(f,zone_corner_southeast zone) -
359 :     sub2(f,zone_corner_northeast zone)) +
360 :     (sub2(f,zone_corner_southwest zone) -
361 :     sub2(f,zone_corner_northwest zone)))
362 :     fun lower_del f =
363 :     0.5 * ((sub2(f,zone_corner_southeast zone) -
364 :     sub2(f,zone_corner_southwest zone)) +
365 :     (sub2(f,zone_corner_northeast zone) -
366 :     sub2(f,zone_corner_northwest zone)))
367 :     val xi = pow(upper_del r',2) + pow(upper_del z',2)
368 :     val eta = pow(lower_del r',2) + pow(lower_del z',2)
369 :     val upper_disc = (upper_del r')*(lower_del w') -
370 :     (upper_del z')*(lower_del u')
371 :     val lower_disc = (upper_del u')*(lower_del z') -
372 :     (upper_del w') * (lower_del r')
373 :     val upper_ubar = if upper_disc<0.0 then upper_disc/xi else 0.0
374 :     val lower_ubar = if lower_disc<0.0 then lower_disc/eta else 0.0
375 :     val gamma = 1.6
376 :     val speed_of_sound = gamma*sub2(p,zone)/sub2(rho',zone)
377 :     val ubar = pow(upper_ubar,2) + pow(lower_ubar,2)
378 :     val viscosity =
379 :     sub2(rho',zone)*(1.5*ubar + 0.5*speed_of_sound*(Math.sqrt ubar))
380 :     val length = Math.sqrt(pow(upper_del r',2) + pow(lower_del r',2))
381 :     val courant_delta = 0.5* sub2(alpha',zone)/(speed_of_sound*length)
382 :     in (viscosity, courant_delta)
383 :     end
384 :     val q' = array2(dimension_all_zones, 0.0)
385 :     val d = array2(dimension_all_zones, 0.0)
386 :     fun reflect_viscosity_cdelta (direction, zone) =
387 :     sub2(q',direction zone) * sub1(qb, sub2(nbc,zone))
388 :     fun do_zones (dir,zone) =
389 :     update2(q',zone,reflect_viscosity_cdelta (dir,zone))
390 :     in
391 :     for_interior_zones (fn zone => let val (qv,dv) = interior_viscosity zone
392 :     in update2(q',zone,qv);
393 :     update2(d,zone,dv)
394 :     end);
395 :     for_south_zones (fn zone => do_zones(north,zone));
396 :     for_east_zones (fn zone => do_zones(west,zone));
397 :     for_west_zones (fn zone => do_zones(east,zone));
398 :     for_north_zones (fn zone => do_zones(south,zone));
399 :     (q', d)
400 :     end
401 :    
402 :     (*
403 :     * Pressure and Energy Polynomial (page 12)
404 :     *)
405 :    
406 :     fun polynomial(G,degree,rho_table,theta_table,rho_value,theta_value) =
407 :     let fun table_search (table, value) =
408 :     let val (low, high) = bounds1 table
409 :     fun search_down i = if value > sub1(table,i-1) then i
410 :     else search_down (i-1)
411 :     in
412 :     if value>sub1(table,high) then high+1
413 :     else if value <= sub1(table,low) then low
414 :     else search_down high
415 :     end
416 :     val rho_index = table_search(rho_table, rho_value)
417 :     val theta_index = table_search(theta_table, theta_value)
418 :     val A = sub2(G, (rho_index, theta_index))
419 :     fun from(n,m) = if n>m then [] else n::from(n+1,m)
420 :     fun f(i,j) = sub2(A,(i,j))*pow(rho_value,i)*pow(theta_value,j)
421 :     in
422 :     sum_list (map (fn i => sum_list(map (fn j => f (i,j)) (from(0,degree))))
423 :     (from (0,degree)))
424 :     end
425 :     fun zonal_pressure (rho_value:real, theta_value:real) =
426 :     let val (G,degree,rho_table,theta_table) =
427 :     extract_pressure_tables_from_constants
428 :     in polynomial(G, degree, rho_table, theta_table, rho_value, theta_value)
429 :     end
430 :    
431 :    
432 :     fun zonal_energy (rho_value, theta_value) =
433 :     let val (G, degree, rho_table, theta_table) =
434 :     extract_energy_tables_from_constants
435 :     in polynomial(G, degree, rho_table, theta_table, rho_value, theta_value)
436 :     end
437 :     val dx = 0.000001
438 :     val tiny = 0.000001
439 :    
440 :    
441 :     fun newton_raphson (f,x) =
442 :     let fun iter (x,fx) =
443 :     if fx > tiny then
444 :     let val fxdx = f(x+dx)
445 :     val denom = fxdx - fx
446 :     in if denom < tiny then iter(x,tiny)
447 :     else iter(x-fx*dx/denom, fxdx)
448 :     end
449 :     else x
450 :     in iter(x, f x)
451 :     end
452 :    
453 :     (*
454 :     * Temperature (page 13-14)
455 :     *)
456 :    
457 :     fun make_temperature(p,epsilon,rho,theta,rho_prime,q_prime) =
458 :     let fun interior_temperature zone =
459 :     let val qkl = sub2(q_prime,zone)
460 :     val rho_kl = sub2(rho,zone)
461 :     val rho_prime_kl = sub2(rho_prime,zone)
462 :     val tau_kl = (1.0 /rho_prime_kl - 1.0/rho_kl)
463 :     fun energy_equation epsilon_kl theta_kl =
464 :     epsilon_kl - zonal_energy(rho_kl,theta_kl)
465 :     val epsilon_0 = sub2(epsilon,zone)
466 :     fun revised_energy pkl = epsilon_0 - (pkl + qkl) * tau_kl
467 :     fun revised_temperature epsilon_kl theta_kl =
468 :     newton_raphson ((energy_equation epsilon_kl), theta_kl)
469 :     fun revised_pressure theta_kl = zonal_pressure(rho_kl, theta_kl)
470 :     val p_0 = sub2(p,zone)
471 :     val theta_0 = sub2(theta,zone)
472 :     val epsilon_1 = revised_energy p_0
473 :     val theta_1 = revised_temperature epsilon_1 theta_0
474 :     val p_1 = revised_pressure theta_1
475 :     val epsilon_2 = revised_energy p_1
476 :     val theta_2 = revised_temperature epsilon_2 theta_1
477 :     in theta_2
478 :     end
479 :     val M = array2(dimension_all_zones, constant_heat_source)
480 :     in
481 :     for_interior_zones
482 :     (fn zone => update2(M, zone, interior_temperature zone));
483 :     M
484 :     end
485 :    
486 :    
487 :     (*
488 :     * Heat conduction
489 :     *)
490 :    
491 :     fun make_cc(alpha_prime, theta_hat) =
492 :     let fun interior_cc zone =
493 :     (0.0001 * pow(sub2(theta_hat,zone),2) *
494 :     (Math.sqrt (abs(sub2(theta_hat,zone)))) / sub2(alpha_prime,zone))
495 :     handle Sqrt => (print (Real.toString (sub2(theta_hat, zone)));
496 :     print ("\nzone =(" ^ Int.toString (#1 zone) ^ "," ^
497 :     Int.toString (#2 zone) ^ ")\n");
498 :     printarray2 theta_hat;
499 :     raise Sqrt)
500 :     val cc = array2(dimension_all_zones, 0.0)
501 :     in
502 :     for_interior_zones(fn zone => update2(cc,zone, interior_cc zone));
503 :     for_south_zones(fn zone => update2(cc,zone, reflect_north zone cc));
504 :     for_west_zones(fn zone => update2(cc,zone,reflect_east zone cc));
505 :     for_east_zones(fn zone => update2(cc,zone,reflect_west zone cc));
506 :     for_north_zones(fn zone => update2(cc,zone, reflect_south zone cc));
507 :     cc
508 :     end
509 :    
510 :     fun make_sigma(deltat, rho_prime, alpha_prime) =
511 :     let fun interior_sigma zone =
512 :     sub2(rho_prime,zone)*sub2(alpha_prime,zone)*specific_heat/ deltat
513 :     val M = array2(dimension_interior_zones, 0.0)
514 :     fun ohandle zone =
515 :     (print (Real.toString (sub2(rho_prime, zone)) ^ " ");
516 :     print (Real.toString (sub2(alpha_prime, zone)) ^ " ");
517 :     print (Real.toString specific_heat ^ " ");
518 :     print (Real.toString deltat ^ "\n");
519 :     raise Overflow)
520 :    
521 :     in if !Control.trace
522 :     then print ("\t\tmake_sigma:deltat = " ^ Real.toString deltat ^ "\n")
523 :     else ();
524 :     (*** for_interior_zones(fn zone => update2(M,zone, interior_sigma zone)) **)
525 :     for_interior_zones(fn zone => (update2(M,zone, interior_sigma zone)
526 :     handle Overflow => ohandle zone));
527 :     M
528 :     end
529 :    
530 :     fun make_gamma ((r_prime,z_prime), cc, succeeding, adjacent) =
531 :     let fun interior_gamma zone =
532 :     let val r1 = sub2(r_prime, zone_corner_southeast zone)
533 :     val z1 = sub2(z_prime, zone_corner_southeast zone)
534 :     val r2 = sub2(r_prime, zone_corner_southeast (adjacent zone))
535 :     val z2 = sub2(z_prime, zone_corner_southeast (adjacent zone))
536 :     val cross_section = 0.5*(r1+r2)*(pow(r1 - r2,2)+pow(z1 - z2,2))
537 :     val (c1,c2) = (sub2(cc, zone), sub2(cc, succeeding zone))
538 :     val specific_conductivity = 2.0 * c1 * c2 / (c1 + c2)
539 :     in cross_section * specific_conductivity
540 :     end
541 :     val M = array2(dimension_all_zones, 0.0)
542 :     in
543 :     for_interior_zones(fn zone => update2(M,zone,interior_gamma zone));
544 :     M
545 :     end
546 :    
547 :     fun make_ab(theta, sigma, Gamma, preceding) =
548 :     let val a = array2(dimension_all_zones, 0.0)
549 :     val b = array2(dimension_all_zones, 0.0)
550 :     fun interior_ab zone =
551 :     let val denom = sub2(sigma, zone) + sub2(Gamma, zone) +
552 :     sub2(Gamma, preceding zone) *
553 :     (1.0 - sub2(a, preceding zone))
554 :     val nume1 = sub2(Gamma,zone)
555 :     val nume2 = sub2(Gamma,preceding zone)*sub2(b,preceding zone) +
556 :     sub2(sigma,zone) * sub2(theta,zone)
557 :     in (nume1/denom, nume2 / denom)
558 :     end
559 :     val f = fn zone => update2(b,zone,sub2(theta,zone))
560 :     in
561 :     for_north_zones f;
562 :     for_south_zones f;
563 :     for_west_zones f;
564 :     for_east_zones f;
565 :     for_interior_zones(fn zone => let val ab = interior_ab zone
566 :     in update2(a,zone,#1 ab);
567 :     update2(b,zone,#2 ab)
568 :     end);
569 :     (a,b)
570 :     end
571 :    
572 :     fun make_theta (a, b, succeeding, int_zones) =
573 :     let val theta = array2(dimension_all_zones, constant_heat_source)
574 :     fun interior_theta zone =
575 :     sub2(a,zone) * sub2(theta,succeeding zone)+ sub2(b,zone)
576 :     in
577 :     int_zones (fn (k,l) => update2(theta, (k,l), interior_theta (k,l)));
578 :     theta
579 :     end
580 :    
581 :     fun compute_heat_conduction(theta_hat, deltat, x', alpha', rho') =
582 :     let val sigma = make_sigma(deltat, rho', alpha')
583 :     val _ = if !Control.trace then print "\tdone make_sigma\n" else ()
584 :    
585 :     val cc = make_cc(alpha', theta_hat)
586 :     val _ = if !Control.trace then print "\tdone make_cc\n" else ()
587 :    
588 :     val Gamma_k = make_gamma( x', cc, north, east)
589 :     val _ = if !Control.trace then print "\tdone make_gamma\n" else ()
590 :    
591 :     val (a_k,b_k) = make_ab(theta_hat, sigma, Gamma_k, north)
592 :     val _ = if !Control.trace then print "\tdone make_ab\n" else ()
593 :    
594 :     val theta_k = make_theta(a_k,b_k,south,for_north_ward_interior_zones)
595 :     val _ = if !Control.trace then print "\tdone make_theta\n" else ()
596 :    
597 :     val Gamma_l = make_gamma(x', cc, west, south)
598 :     val _ = if !Control.trace then print "\tdone make_gamma\n" else ()
599 :    
600 :     val (a_l,b_l) = make_ab(theta_k, sigma, Gamma_l, west)
601 :     val _ = if !Control.trace then print "\tdone make_ab\n" else ()
602 :    
603 :     val theta_l = make_theta(a_l,b_l,east,for_west_ward_interior_zones)
604 :     val _ = if !Control.trace then print "\tdone make_theta\n" else ()
605 :     in (theta_l, Gamma_k, Gamma_l)
606 :     end
607 :    
608 :    
609 :     (*
610 :     * Final Pressure and Energy calculation
611 :     *)
612 :     fun make_pressure(rho', theta') =
613 :     let val p = array2(dimension_all_zones, 0.0)
614 :     fun boundary_p(direction, zone) =
615 :     sub1(pbb, sub2(nbc, zone)) +
616 :     sub1(pb,sub2(nbc,zone)) * sub2(p, direction zone)
617 :     in
618 :     for_interior_zones
619 :     (fn zone => update2(p,zone,zonal_pressure(sub2(rho',zone),
620 :     sub2(theta',zone))));
621 :     for_south_zones(fn zone => update2(p,zone,boundary_p(north,zone)));
622 :     for_east_zones(fn zone => update2(p,zone,boundary_p(west,zone)));
623 :     for_west_zones(fn zone => update2(p,zone,boundary_p(east,zone)));
624 :     for_north_zones(fn zone => update2(p,zone,boundary_p(south,zone)));
625 :     p
626 :     end
627 :    
628 :     fun make_energy(rho', theta') =
629 :     let val epsilon' = array2(dimension_all_zones, 0.0)
630 :     in
631 :     for_interior_zones
632 :     (fn zone => update2(epsilon', zone, zonal_energy(sub2(rho',zone),
633 :     sub2(theta',zone))));
634 :     for_south_zones
635 :     (fn zone => update2(epsilon',zone, reflect_north zone epsilon'));
636 :     for_west_zones
637 :     (fn zone => update2(epsilon',zone, reflect_east zone epsilon'));
638 :     for_east_zones
639 :     (fn zone => update2(epsilon',zone, reflect_west zone epsilon'));
640 :     for_north_zones
641 :     (fn zone => update2(epsilon',zone, reflect_south zone epsilon'));
642 :     epsilon'
643 :     end
644 :    
645 :    
646 :     (*
647 :     * Energy Error Calculation (page 20)
648 :     *)
649 :    
650 :     fun compute_energy_error ((u',w'),(r',z'),p',q',epsilon',theta',rho',alpha',
651 :     Gamma_k,Gamma_l,deltat) =
652 :     let fun mass zone = sub2(rho',zone) * sub2(alpha',zone):real
653 :     val internal_energy =
654 :     sum_list (map_interior_zones (fn z => sub2(epsilon',z)*(mass z)))
655 :     fun kinetic node =
656 :     let val average_mass = 0.25*((mass (zone_A node)) +
657 :     (mass (zone_B node)) +
658 :     (mass (zone_C node)) +
659 :     (mass (zone_D node)))
660 :     val v_square = pow(sub2(u',node),2) + pow(sub2(w',node),2)
661 :     in 0.5 * average_mass * v_square
662 :     end
663 :     val kinetic_energy = sum_list (map_interior_nodes kinetic)
664 :     fun work_done (node1, node2) =
665 :     let val (r1, r2) = (sub2(r',node1), sub2(r',node2))
666 :     val (z1, z2) = (sub2(z',node1), sub2(z',node2))
667 :     val (u1, u2) = (sub2(p',node1), sub2(p',node2))
668 :     val (w1, w2) = (sub2(z',node1), sub2(z',node2))
669 :     val (p1, p2) = (sub2(p',node1), sub2(p',node2))
670 :     val (q1, q2) = (sub2(q',node1), sub2(q',node2))
671 :     val force = 0.5*(p1+p2+q1+q2)
672 :     val radius = 0.5* (r1+r2)
673 :     val area = 0.5* ((r1-r2)*(u1-u2) - (z1-z2)*(w1-w2))
674 :     in force * radius * area * deltat
675 :     end
676 :    
677 :     fun from(n,m) = if n > m then [] else n::from(n+1,m)
678 :     val north_line =
679 :     map (fn l => (west(kmin,l),(kmin,l))) (from(lmin+1,lmax))
680 :     val south_line =
681 :     map (fn l => (west(kmax,l),(kmax,l))) (from(lmin+1,lmax))
682 :     val east_line =
683 :     map (fn k => (south(k,lmax),(k,lmax))) (from(kmin+1,kmax))
684 :     val west_line =
685 :     map (fn k => (south(k,lmin+1),(k,lmin+1))) (from(kmin+1,kmax))
686 :    
687 :     val w1 = sum_list (map work_done north_line)
688 :     val w2 = sum_list (map work_done south_line)
689 :     val w3 = sum_list (map work_done east_line)
690 :     val w4 = sum_list (map work_done west_line)
691 :     val boundary_work = w1 + w2 + w3 + w4
692 :    
693 :     fun heat_flow Gamma (zone1,zone2) =
694 :     deltat * sub2(Gamma, zone1) * (sub2(theta',zone1) - sub2(theta',zone2))
695 :    
696 :     val north_flow =
697 :     let val k = kmin+1
698 :     in map (fn l => (north(k,l),(k,l))) (from(lmin+1,lmax))
699 :     end
700 :     val south_flow =
701 :     let val k = kmax
702 :     in map (fn l => (south(k,l),(k,l))) (from(lmin+2,lmax-1))
703 :     end
704 :     val east_flow =
705 :     let val l = lmax
706 :     in map (fn k => (east(k,l),(k,l))) (from(kmin+2,kmax))
707 :     end
708 :     val west_flow =
709 :     let val l = lmin+1
710 :     in map (fn k => (west(k,l),(k,l))) (from(kmin+2,kmax))
711 :     end
712 :    
713 :     val h1 = sum_list (map (heat_flow Gamma_k) north_flow)
714 :     val h2 = sum_list (map (heat_flow Gamma_k) south_flow)
715 :     val h3 = sum_list (map (heat_flow Gamma_l) east_flow)
716 :     val h4 = sum_list (map (heat_flow Gamma_l) west_flow)
717 :     val boundary_heat = h1 + h2 + h3 + h4
718 :     in
719 :     internal_energy + kinetic_energy - boundary_heat - boundary_work
720 :     end
721 :    
722 :     fun compute_time_step(d, theta_hat, theta') =
723 :     let val deltat_courant =
724 :     min_list (map_interior_zones (fn zone => sub2(d,zone)))
725 :     val deltat_conduct =
726 :     max_list (map_interior_zones
727 :     (fn z => (abs(sub2(theta_hat,z) - sub2(theta', z))/
728 :     sub2(theta_hat,z))))
729 :     val deltat_minimum = min (deltat_courant, deltat_conduct)
730 :     in min (deltat_maximum, deltat_minimum)
731 :     end
732 :    
733 :    
734 :     fun compute_initial_state () =
735 :     let
736 :     val v = (all_zero_nodes, all_zero_nodes)
737 :     val x = let fun interior_position (k,l) =
738 :     let val pi = 3.1415926535898
739 :     val rp = real (lmax - lmin)
740 :     val z1 = real(10 + k - kmin)
741 :     val zz = (~0.5 + real(l - lmin) / rp) * pi
742 :     in (z1 * Math.cos zz, z1 * Math.sin zz)
743 :     end
744 :     in make_position_matrix interior_position
745 :     end
746 :     val (alpha,s) =
747 :     let val (alpha_prime,s_prime) =
748 :     let val A = array2(dimension_all_zones, 0.0)
749 :     val S = array2(dimension_all_zones, 0.0)
750 :     fun reflect_area_vol f = (f A, f S)
751 :    
752 :     fun u2 (f,z) =
753 :     let val (a,s) = reflect_area_vol(f z)
754 :     in update2(A,z,a);
755 :     update2(S,z,s)
756 :     end
757 :     in
758 :     for_interior_zones
759 :     (fn z => let val (a,s) = zone_area_vol(x, z)
760 :     in update2(A,z,a);
761 :     update2(S,z,s)
762 :     end);
763 :     for_south_zones (fn z => u2 (reflect_north, z));
764 :     for_east_zones (fn z => u2 (reflect_west, z));
765 :     for_west_zones (fn z => u2 (reflect_east, z));
766 :     for_north_zones (fn z => u2 (reflect_south, z));
767 :     (A,S)
768 :     end
769 :     in (alpha_prime,s_prime)
770 :     end
771 :     val rho = let val R = array2(dimension_all_zones, 0.0)
772 :     in for_all_zones (fn z => update2(R,z,1.4)); R
773 :     end
774 :     val theta =
775 :     let val T = array2(dimension_all_zones, constant_heat_source)
776 :     in for_interior_zones(fn z => update2(T,z,0.0001));
777 :     T
778 :     end
779 :     val p = make_pressure(rho, theta)
780 :     val q = all_zero_zones
781 :     val epsilon = make_energy(rho, theta)
782 :     val deltat = 0.01
783 :     val c = 0.0
784 :     in
785 :     (v,x,alpha,s,rho,p,q,epsilon,theta,deltat,c)
786 :     end
787 :    
788 :    
789 :     fun compute_next_state state =
790 :     let
791 :     val (v,x,alpha,s,rho,p,q,epsilon,theta,deltat,c) = state
792 :     val v' = make_velocity (v, x, p, q, alpha, rho, deltat)
793 :     val _ = if !Control.trace then print "done make_velocity\n" else ()
794 :    
795 :     val x' = make_position(x,deltat,v')
796 :     handle Overflow =>(printarray2 (#1 v');
797 :     printarray2 (#2 v');
798 :     raise Overflow)
799 :     val _ = if !Control.trace then print "done make_position\n" else ()
800 :    
801 :     val (alpha',rho',s') = make_area_density_volume (rho, s , x')
802 :     val _ = if !Control.trace then print "done make_area_density_volume\n"
803 :     else ()
804 :    
805 :     val (q',d) = make_viscosity (p, v', x', alpha', rho')
806 :     val _ = if !Control.trace then print "done make_viscosity\n" else ()
807 :    
808 :     val theta_hat = make_temperature (p, epsilon, rho, theta, rho', q')
809 :     val _ = if !Control.trace then print "done make_temperature\n" else ()
810 :    
811 :     val (theta',Gamma_k,Gamma_l) =
812 :     compute_heat_conduction (theta_hat, deltat, x', alpha', rho')
813 :     val _ = if !Control.trace then print "done compute_heat_conduction\n"
814 :     else ()
815 :    
816 :     val p' = make_pressure(rho', theta')
817 :     val _ = if !Control.trace then print "done make_pressure\n" else ()
818 :    
819 :     val epsilon' = make_energy (rho', theta')
820 :     val _ = if !Control.trace then print "done make_energy\n" else ()
821 :    
822 :     val c' = compute_energy_error (v', x', p', q', epsilon', theta', rho',
823 :     alpha', Gamma_k, Gamma_l, deltat)
824 :     val _ = if !Control.trace then print "done compute_energy_error\n"
825 :     else ()
826 :    
827 :     val deltat' = compute_time_step (d, theta_hat, theta')
828 :     val _ = if !Control.trace then print "done compute_time_step\n\n" else ()
829 :     in
830 :     (v',x',alpha',s',rho',p',q', epsilon',theta',deltat',c')
831 :     end
832 :    
833 :     fun runit () =
834 :     let fun iter (i,state) = if i = 0 then state
835 :     else (print ".";
836 :     iter(i-1, compute_next_state state))
837 :     in iter(step_count, compute_initial_state())
838 :     end
839 :    
840 :     fun print_state ((v1,v2),(r,z),alpha,s,rho,p,q,epsilon,theta,deltat,c) = (
841 :     print "Velocity matrices = \n";
842 :     printarray2 v1; print "\n\n";
843 :     printarray2 v2;
844 :    
845 :     print "\n\nPosition matrices = \n";
846 :     printarray2 r; print "\n\n";
847 :     printarray2 z;
848 :    
849 :     print "\n\nalpha = \n";
850 :     printarray2 alpha;
851 :    
852 :     print "\n\ns = \n";
853 :     printarray2 s;
854 :    
855 :     print "\n\nrho = \n";
856 :     printarray2 rho;
857 :    
858 :     print "\n\nPressure = \n";
859 :     printarray2 p;
860 :    
861 :     print "\n\nq = \n";
862 :     printarray2 q;
863 :    
864 :     print "\n\nepsilon = \n";
865 :     printarray2 epsilon;
866 :    
867 :     print "\n\ntheta = \n";
868 :     printarray2 theta;
869 :    
870 :     print ("delatat = " ^ Real.toString (deltat : real)^ "\n");
871 :     print ("c = " ^ Real.toString (c : real) ^ "\n"))
872 :    
873 :     fun testit outstrm = print_state (runit())
874 :    
875 :     fun doit () = let
876 :     val (_, _, _, _, _, _, _, _, _, delta', c') = runit()
877 :     val delta = Real.trunc delta'
878 :     val c = Real.trunc (c' * 10000.0)
879 :     in
880 :     if (c = 6787 andalso delta = ~33093)
881 :     then ()
882 :     else TextIO.output (TextIO.stdErr, "*** ERROR ***\n")
883 :     end
884 :    
885 :     end; (* functor Simple *)

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