1 (* From the SML
/NJ benchmark suite
. *)
11 val array
: (int*int) * 'a
-> 'a array2
12 val sub
: 'a array2
* (int*int) -> 'a
13 val update
: 'a array2
* (int*int) * 'a
-> unit
14 val length
: 'a array2
-> (int*int)
18 type 'a array2
= {size
: (int*int), value
: 'a Array
.array
}
19 exception Subscript
= Subscript
20 fun index ((i1
:int,i2
:int),(s1
,s2
)) =
21 if i1
>=0 andalso i1
<s1
andalso i2
>=0 andalso i2
<s2
then i1
*s2
+i2
23 fun array(bnds
as (i1
,i2
), v
) = {size
=bnds
, value
=Array
.array(i1
*i2
, v
)}
24 fun op sub ({size
,value
}, indx
) = Array
.sub(value
, index(indx
,size
))
25 fun update ({size
=size
,value
=A
},i
,v
) = Array
.update(A
,index(i
,size
),v
)
26 fun length
{size
=size
,value
=A
} = size
31 val doit
: int -> unit
32 val testit
: TextIO.outstream
-> unit
37 functor Simple(val grid_max
: int val step_count
: int) : BMARK
=
40 fun fold f
[] = (fn b
=> b
)
41 | fold
f (a
::r
) = (fn b
=> let fun f2(e
,[]) = f(e
,b
)
42 |
f2(e
,a
::r
) = f(e
,f2(a
,r
))
47 fun min (x
:real,y
:real) = if x
<y
then x
else y
48 fun max (x
:real,y
:real) = if x
<y
then y
else x
52 fun max_list
[] = raise MaxList | max_list l
= fold max
l (hd l
)
53 fun min_list
[] = raise MinList | min_list l
= fold min
l (hd l
)
54 fun sum_list
[] = raise SumList
55 |
sum_list (l
:real list
) = fold (op +) l
0.0
57 fun for
{from
=start
:int,step
=delta
:int, to
=endd
:int} body
=
58 if delta
>0 andalso endd
>=start
then
59 let fun f x
= if x
> endd
then () else (body x
; f(x
+delta
))
62 else if endd
<=start
then
63 let fun f x
= if x
< endd
then () else (body x
; f(x
+delta
))
67 fun from(n
,m
) = if n
>m
then [] else n
::from(n
+1,m
)
69 |
flatten (x
::xs
) = x @ flatten xs
70 fun pow(x
:real,y
:int) = if y
= 0 then 1.0 else x
* pow(x
,y
-1)
71 fun array2(bounds
as ((l1
,u1
),(l2
,u2
)),v
) =
72 (Array2
.array((u1
-l1
+1, u2
-l2
+1),v
), bounds
)
73 fun sub2((A
,((lb1
:int,ub1
:int),(lb2
:int,ub2
:int))),(k
,l
)) =
74 Array2
.sub(A
, (k
-lb1
, l
-lb2
))
75 fun update2((A
,((lb1
,_
),(lb2
,_
))),(k
,l
), v
) = Array2
.update(A
,(k
-lb1
,l
-lb2
),v
)
77 fun printarray2 (A
as (M
:real Array2
.array2
,((l1
,u1
),(l2
,u2
)))) =
78 for
{from
=l1
,step
=1,to
=u1
} (fn i
=>
80 for
{from
=l2
,step
=1,to
=u2
-1} (fn j
=>
81 print (Real.toString (sub2(A
,(i
,j
))) ^
", "));
82 print (Real.toString (sub2(A
,(i
,u2
))) ^
"]\n")))
83 fun array1((l
,u
),v
) = (Array
.array(u
-l
+1,v
),(l
,u
))
84 fun sub1((A
,(l
:int,u
:int)),i
:int) = Array
.sub(A
,i
-l
)
85 fun update1((A
,(l
,_
)),i
,v
) = Array
.update(A
,i
-l
,v
)
89 * Specification
of the state variable computation
91 val grid_size
= ((2,grid_max
), (2,grid_max
))
93 fun north (k
,l
) = (k
-1,l
)
94 fun south (k
,l
) = (k
+1,l
)
96 fun east (k
,l
) = (k
,l
+1)
97 fun west (k
,l
) = (k
,l
-1)
99 val northeast
= north
o east
100 val southeast
= south
o east
101 val northwest
= north
o west
102 val southwest
= south
o west
104 fun farnorth x
= (north
o north
) x
105 fun farsouth x
= (south
o south
) x
106 fun fareast x
= (east
o east
) x
107 fun farwest x
= (west
o west
) x
109 fun zone_A(k
,l
) = (k
,l
)
110 fun zone_B(k
,l
) = (k
+1,l
)
112 fun zone_C(k
,l
) = (k
+1,l
+1)
113 fun zone_D(k
,l
) = (k
,l
+1)
115 val zone_corner_northeast
= north
116 val zone_corner_northwest
= northwest
117 fun zone_corner_southeast zone
= zone
118 val zone_corner_southwest
= west
120 val ((kmin
,kmax
),(lmin
,lmax
)) = grid_size
121 val dimension_all_nodes
= ((kmin
-1,kmax
+1),(lmin
-1,lmax
+1))
122 fun for_all_nodes f
=
123 for
{from
=kmin
-1, step
=1, to
=kmax
+1} (fn k
=>
124 for
{from
=lmin
-1, step
=1, to
=lmax
+1} (fn l
=> f k l
))
126 val dimension_interior_nodes
= ((kmin
,kmax
),(lmin
,lmax
))
127 fun for_interior_nodes f
=
128 for
{from
=kmin
, step
=1, to
=kmax
} (fn k
=>
129 for
{from
=lmin
, step
=1, to
=lmax
} (fn l
=> f k l
))
131 val dimension_all_zones
= ((kmin
,kmax
+1),(lmin
,lmax
+1))
132 fun for_all_zones f
=
133 for
{from
=kmin
, step
=1, to
=kmax
+1} (fn k
=>
134 for
{from
=lmin
, step
=1, to
=lmax
+1} (fn l
=> f (k
,l
)))
136 val dimension_interior_zones
= ((kmin
+1,kmax
),(lmin
+1,lmax
))
137 fun for_interior_zones f
=
138 for
{from
=kmin
+1, step
=1, to
=kmax
} (fn k
=>
139 for
{from
=lmin
+1, step
=1, to
=lmax
} (fn l
=> f (k
,l
)))
141 fun map_interior_nodes f
=
142 flatten(map (fn k
=> (map (fn l
=> f (k
,l
))
145 fun map_interior_zones f
=
146 flatten(map (fn k
=> (map (fn l
=> f (k
,l
))
147 (from(lmin
+1,lmax
))))
150 fun for_north_ward_interior_zones f
=
151 for
{from
=kmax
, step
= ~
1, to
=kmin
+1} (fn k
=>
152 for
{from
=lmin
+1, step
=1, to
=lmax
} (fn l
=> f (k
,l
)))
153 fun for_west_ward_interior_zones f
=
154 for
{from
=kmin
+1, step
=1, to
=kmax
} (fn k
=>
155 for
{from
=lmax
, step
= ~
1, to
=lmin
+1} (fn l
=> f (k
,l
)))
158 fun for_north_zones f
= for
{from
=lmin
, step
=1, to
=lmax
+1} (fn l
=> f (kmin
,l
))
159 fun for_south_zones f
= for
{from
=lmin
+1, step
=1, to
=lmax
} (fn l
=> f (kmax
+1,l
))
160 fun for_east_zones f
= for
{from
=kmin
+1, step
=1, to
=kmax
+1}(fn k
=> f (k
,lmax
+1))
161 fun for_west_zones f
= for
{from
=kmin
+1, step
=1, to
=kmax
+1}(fn k
=> f (k
,lmin
))
163 fun reflect dir node A
= sub2(A
, dir node
)
164 val reflect_north
= fn x
=> reflect north x
165 val reflect_south
= fn x
=> reflect south x
166 val reflect_east
= fn x
=> reflect east x
167 val reflect_west
= fn x
=> reflect west x
169 fun for_north_nodes f
=
170 for
{from
=lmin
, step
=1, to
=lmax
-1} (fn l
=> f (kmin
-1,l
))
171 fun for_south_nodes f
=
172 for
{from
=lmin
, step
=1, to
=lmax
-1} (fn l
=> f (kmax
+1,l
))
173 fun for_east_nodes f
=
174 for
{from
=kmin
, step
=1, to
=kmax
-1} (fn k
=> f (k
,lmax
+1))
175 fun for_west_nodes f
=
176 for
{from
=kmin
, step
=1, to
=kmax
-1} (fn k
=> f (k
,lmin
-1))
178 val north_east_corner
= (kmin
-1,lmax
+1)
179 val north_west_corner
= (kmin
-1,lmin
-1)
180 val south_east_corner
= (kmax
+1,lmax
+1)
181 val south_west_corner
= (kmax
+1,lmin
-1)
183 val west_of_north_east
= (kmin
-1, lmax
)
184 val west_of_south_east
= (kmax
+1, lmax
)
185 val north_of_south_east
= (kmax
, lmax
+1)
186 val north_of_south_west
= (kmax
, lmin
-1)
191 * Initialization
of parameters
193 val constant_heat_source
= 0.0
194 val deltat_maximum
= 0.01
195 val specific_heat
= 0.1
196 val p_coeffs
= let val M
= array2(((0,2),(0,2)), 0.0)
197 in update2(M
, (1,1), 0.06698); M
199 val e_coeffs
= let val M
= array2(((0,2),(0,2)), 0.0)
200 in update2(M
, (0,1), 0.1); M
202 val p_poly
= array2(((1,4),(1,5)),p_coeffs
)
204 val e_poly
= array2(((1,4),(1,5)), e_coeffs
)
206 val rho_table
= let val V
= array1((1,3), 0.0)
211 val theta_table
= let val V
= array1((1,4), 0.0)
218 val extract_energy_tables_from_constants
= (e_poly
,2,rho_table
,theta_table
)
219 val extract_pressure_tables_from_constants
= (p_poly
,2,rho_table
,theta_table
)
221 val nbc
= let val M
= array2(dimension_all_zones
, 1)
222 in for
{from
=lmin
+1,step
=1,to
=lmax
} (fn j
=> update2(M
,(kmax
+1, j
),2));
223 update2(M
,(kmin
,lmin
),4);
224 update2(M
,(kmin
,lmax
+1),4);
225 update2(M
,(kmax
+1,lmin
),4);
226 update2(M
,(kmax
+1,lmax
+1),4);
229 val pbb
= let val A
= array1((1,4), 0.0)
230 in update1(A
,2,6.0); A
232 val pb
= let val A
= array1((1,4), 1.0)
233 in update1(A
,2,0.0); update1(A
,3,0.0); A
237 val all_zero_nodes
= array2(dimension_all_nodes
, 0.0)
239 val all_zero_zones
= array2(dimension_all_zones
, 0.0)
243 * Positional Coordinates
. (page
9-10)
246 fun make_position_matrix interior_function
=
247 let val r
' = array2(dimension_all_nodes
, 0.0)
248 val z
' = array2(dimension_all_nodes
, 0.0)
249 fun boundary_position (rx
,zx
,ry
,zy
,ra
,za
) =
250 let val (rax
, zax
) = (ra
- rx
, za
- zx
)
251 val (ryx
, zyx
) = (ry
- rx
, zy
- zx
)
252 val omega
= 2.0*(rax
*ryx
+ zax
*zyx
)/(ryx
*ryx
+ zyx
*zyx
)
253 val rb
= rx
- rax
+ omega
*ryx
254 val zb
= zx
- zax
+ omega
*zyx
258 fun reflect_node (x_dir
, y_dir
, a_dir
, node
) =
259 let val rx
= reflect x_dir node r
'
260 val zx
= reflect x_dir node z
'
261 val ry
= reflect y_dir node r
'
262 val zy
= reflect y_dir node z
'
263 val ra
= reflect a_dir node r
'
264 val za
= reflect a_dir node z
'
265 in boundary_position (rx
, zx
, ry
, zy
, ra
, za
)
267 fun u2 (rv
,zv
) n
= (update2(r
',n
,rv
); update2(z
',n
,zv
))
269 for_interior_nodes (fn k
=> fn l
=> u2 (interior_function (k
,l
)) (k
,l
));
270 for_north_nodes(fn n
=> u2 (reflect_node(south
,southeast
,farsouth
,n
)) n
);
271 for_south_nodes (fn n
=> u2(reflect_node(north
,northeast
,farnorth
,n
)) n
);
272 for_east_nodes (fn n
=> u2(reflect_node(west
, southwest
, farwest
, n
)) n
);
273 for_west_nodes (fn n
=> u2(reflect_node(east
, southeast
, fareast
, n
)) n
);
274 u2 (reflect_node(south
, southwest
, farsouth
, west_of_north_east
))
276 u2 (reflect_node(north
, northwest
, farnorth
, west_of_south_east
))
278 u2 (reflect_node(west
, northwest
, farwest
, north_of_south_east
))
280 u2 (reflect_node(east
, northeast
, fareast
, north_of_south_west
))
282 u2 (reflect_node(southwest
, west
, farwest
, north_east_corner
))
284 u2 (reflect_node(northwest
, west
, farwest
, south_east_corner
))
286 u2 (reflect_node(southeast
, south
, farsouth
, north_west_corner
))
288 u2 (reflect_node(northeast
, east
, fareast
, south_west_corner
))
296 * Physical Properties
of a
Zone (page
10)
298 fun zone_area_vol ((r
,z
), zone
) =
299 let val (r1
,z1
)=(sub2(r
,zone_corner_southwest zone
),
300 sub2(z
,zone_corner_southwest zone
))
301 val (r2
,z2
)=(sub2(r
,zone_corner_southeast zone
),
302 sub2(z
,zone_corner_southeast zone
))
303 val (r3
,z3
)=(sub2(r
,zone_corner_northeast zone
),
304 sub2(z
,zone_corner_northeast zone
))
305 val (r4
,z4
)=(sub2(r
,zone_corner_northwest zone
),
306 sub2(z
,zone_corner_northwest zone
))
307 val area1
= (r2
-r1
)*(z3
-z1
) - (r3
-r2
)*(z3
-z2
)
308 val radius1
= 0.3333 *(r1
+r2
+r3
)
309 val volume1
= area1
* radius1
310 val area2
= (r3
-r1
)*(z4
-z3
) - (r4
-r3
)*(z3
-z1
)
311 val radius2
= 0.3333 *(r1
+r3
+r4
)
312 val volume2
= area2
* radius2
313 in (area1
+area2
, volume1
+volume2
)
319 fun make_velocity((u
,w
),(r
,z
),p
,q
,alpha
,rho
,delta_t
) =
320 let fun line_integral (p
,z
,node
) : real =
321 sub2(p
,zone_A node
)*(sub2(z
,west node
) - sub2(z
,north node
)) +
322 sub2(p
,zone_B node
)*(sub2(z
,south node
) - sub2(z
,west node
)) +
323 sub2(p
,zone_C node
)*(sub2(z
,east node
) - sub2(z
,south node
)) +
324 sub2(p
,zone_D node
)*(sub2(z
,north node
) - sub2(z
,east node
))
325 fun regional_mass node
=
326 0.5 * (sub2(rho
, zone_A node
)*sub2(alpha
,zone_A node
) +
327 sub2(rho
, zone_B node
)*sub2(alpha
,zone_B node
) +
328 sub2(rho
, zone_C node
)*sub2(alpha
,zone_C node
) +
329 sub2(rho
, zone_D node
)*sub2(alpha
,zone_D node
))
331 let val d
= regional_mass node
332 val n1
= ~
(line_integral(p
,z
,node
)) - line_integral(q
,z
,node
)
333 val n2
= line_integral(p
,r
,node
) + line_integral(q
,r
,node
)
336 in (sub2(u
,node
)+delta_t
*u_dot
, sub2(w
,node
)+delta_t
*w_dot
)
338 val U
= array2(dimension_interior_nodes
,0.0)
339 val W
= array2(dimension_interior_nodes
,0.0)
340 in for_interior_nodes (fn k
=> fn l
=> let val (uv
,wv
) = velocity (k
,l
)
341 in update2(U
,(k
,l
),uv
);
349 fun make_position ((r
,z
),delta_t
,(u
',w
')) =
350 let fun interior_position node
=
351 (sub2(r
,node
) + delta_t
*sub2(u
',node
),
352 sub2(z
,node
) + delta_t
*sub2(w
',node
))
353 in make_position_matrix interior_position
357 fun make_area_density_volume(rho
, s
, x
') =
358 let val alpha
' = array2(dimension_all_zones
, 0.0)
359 val s
' = array2(dimension_all_zones
, 0.0)
360 val rho
' = array2(dimension_all_zones
, 0.0)
361 fun interior_area zone
=
362 let val (area
, vol
) = zone_area_vol (x
', zone
)
363 val density
= sub2(rho
,zone
)*sub2(s
,zone
) / vol
364 in (area
,vol
,density
)
366 fun reflect_area_vol_density reflect_function
=
367 (reflect_function alpha
',reflect_function s
',reflect_function rho
')
368 fun update_asr (zone
,(a
,s
,r
)) = (update2(alpha
',zone
,a
);
370 update2(rho
',zone
,r
))
371 fun r_area_vol_den (reflect_dir
,zone
) =
372 let val asr
= reflect_area_vol_density (reflect_dir zone
)
373 in update_asr(zone
, asr
)
376 for_interior_zones (fn zone
=> update_asr(zone
, interior_area zone
));
377 for_south_zones (fn zone
=> r_area_vol_den(reflect_north
, zone
));
378 for_east_zones (fn zone
=> r_area_vol_den(reflect_west
, zone
));
379 for_west_zones (fn zone
=> r_area_vol_den(reflect_east
, zone
));
380 for_north_zones (fn zone
=> r_area_vol_den(reflect_south
, zone
));
386 * Artifical
Viscosity (page
11)
388 fun make_viscosity(p
,(u
',w
'),(r
',z
'), alpha
',rho
') =
389 let fun interior_viscosity zone
=
390 let fun upper_del f
=
391 0.5 * ((sub2(f
,zone_corner_southeast zone
) -
392 sub2(f
,zone_corner_northeast zone
)) +
393 (sub2(f
,zone_corner_southwest zone
) -
394 sub2(f
,zone_corner_northwest zone
)))
396 0.5 * ((sub2(f
,zone_corner_southeast zone
) -
397 sub2(f
,zone_corner_southwest zone
)) +
398 (sub2(f
,zone_corner_northeast zone
) -
399 sub2(f
,zone_corner_northwest zone
)))
400 val xi
= pow(upper_del r
',2) + pow(upper_del z
',2)
401 val eta
= pow(lower_del r
',2) + pow(lower_del z
',2)
402 val upper_disc
= (upper_del r
')*(lower_del w
') -
403 (upper_del z
')*(lower_del u
')
404 val lower_disc
= (upper_del u
')*(lower_del z
') -
405 (upper_del w
') * (lower_del r
')
406 val upper_ubar
= if upper_disc
<0.0 then upper_disc
/xi
else 0.0
407 val lower_ubar
= if lower_disc
<0.0 then lower_disc
/eta
else 0.0
409 val speed_of_sound
= gamma
*sub2(p
,zone
)/sub2(rho
',zone
)
410 val ubar
= pow(upper_ubar
,2) + pow(lower_ubar
,2)
412 sub2(rho
',zone
)*(1.5*ubar
+ 0.5*speed_of_sound
*(Math
.sqrt ubar
))
413 val length
= Math
.sqrt(pow(upper_del r
',2) + pow(lower_del r
',2))
414 val courant_delta
= 0.5* sub2(alpha
',zone
)/(speed_of_sound
*length
)
415 in (viscosity
, courant_delta
)
417 val q
' = array2(dimension_all_zones
, 0.0)
418 val d
= array2(dimension_all_zones
, 0.0)
419 fun reflect_viscosity_cdelta (direction
, zone
) =
420 sub2(q
',direction zone
) * sub1(qb
, sub2(nbc
,zone
))
421 fun do_zones (dir
,zone
) =
422 update2(q
',zone
,reflect_viscosity_cdelta (dir
,zone
))
424 for_interior_zones (fn zone
=> let val (qv
,dv
) = interior_viscosity zone
425 in update2(q
',zone
,qv
);
428 for_south_zones (fn zone
=> do_zones(north
,zone
));
429 for_east_zones (fn zone
=> do_zones(west
,zone
));
430 for_west_zones (fn zone
=> do_zones(east
,zone
));
431 for_north_zones (fn zone
=> do_zones(south
,zone
));
436 * Pressure
and Energy
Polynomial (page
12)
439 fun polynomial(G
,degree
,rho_table
,theta_table
,rho_value
,theta_value
) =
440 let fun table_search (table
, value
) =
441 let val (low
, high
) = bounds1 table
442 fun search_down i
= if value
> sub1(table
,i
-1) then i
443 else search_down (i
-1)
445 if value
>sub1(table
,high
) then high
+1
446 else if value
<= sub1(table
,low
) then low
447 else search_down high
449 val rho_index
= table_search(rho_table
, rho_value
)
450 val theta_index
= table_search(theta_table
, theta_value
)
451 val A
= sub2(G
, (rho_index
, theta_index
))
452 fun from(n
,m
) = if n
>m
then [] else n
::from(n
+1,m
)
453 fun f(i
,j
) = sub2(A
,(i
,j
))*pow(rho_value
,i
)*pow(theta_value
,j
)
455 sum_list (map (fn i
=> sum_list(map (fn j
=> f (i
,j
)) (from(0,degree
))))
458 fun zonal_pressure (rho_value
:real, theta_value
:real) =
459 let val (G
,degree
,rho_table
,theta_table
) =
460 extract_pressure_tables_from_constants
461 in polynomial(G
, degree
, rho_table
, theta_table
, rho_value
, theta_value
)
465 fun zonal_energy (rho_value
, theta_value
) =
466 let val (G
, degree
, rho_table
, theta_table
) =
467 extract_energy_tables_from_constants
468 in polynomial(G
, degree
, rho_table
, theta_table
, rho_value
, theta_value
)
474 fun newton_raphson (f
,x
) =
475 let fun iter (x
,fx
) =
477 let val fxdx
= f(x
+dx
)
478 val denom
= fxdx
- fx
479 in if denom
< tiny
then iter(x
,tiny
)
480 else iter(x
-fx
*dx
/denom
, fxdx
)
487 * Temperature (page
13-14)
490 fun make_temperature(p
,epsilon
,rho
,theta
,rho_prime
,q_prime
) =
491 let fun interior_temperature zone
=
492 let val qkl
= sub2(q_prime
,zone
)
493 val rho_kl
= sub2(rho
,zone
)
494 val rho_prime_kl
= sub2(rho_prime
,zone
)
495 val tau_kl
= (1.0 /rho_prime_kl
- 1.0/rho_kl
)
496 fun energy_equation epsilon_kl theta_kl
=
497 epsilon_kl
- zonal_energy(rho_kl
,theta_kl
)
498 val epsilon_0
= sub2(epsilon
,zone
)
499 fun revised_energy pkl
= epsilon_0
- (pkl
+ qkl
) * tau_kl
500 fun revised_temperature epsilon_kl theta_kl
=
501 newton_raphson ((energy_equation epsilon_kl
), theta_kl
)
502 fun revised_pressure theta_kl
= zonal_pressure(rho_kl
, theta_kl
)
503 val p_0
= sub2(p
,zone
)
504 val theta_0
= sub2(theta
,zone
)
505 val epsilon_1
= revised_energy p_0
506 val theta_1
= revised_temperature epsilon_1 theta_0
507 val p_1
= revised_pressure theta_1
508 val epsilon_2
= revised_energy p_1
509 val theta_2
= revised_temperature epsilon_2 theta_1
512 val M
= array2(dimension_all_zones
, constant_heat_source
)
515 (fn zone
=> update2(M
, zone
, interior_temperature zone
));
524 fun make_cc(alpha_prime
, theta_hat
) =
525 let fun interior_cc zone
=
526 (0.0001 * pow(sub2(theta_hat
,zone
),2) *
527 (Math
.sqrt (abs(sub2(theta_hat
,zone
)))) / sub2(alpha_prime
,zone
))
528 handle Sqrt
=> (print (Real.toString (sub2(theta_hat
, zone
)));
529 print ("\nzone =(" ^
Int.toString (#
1 zone
) ^
"," ^
530 Int.toString (#
2 zone
) ^
")\n");
531 printarray2 theta_hat
;
533 val cc
= array2(dimension_all_zones
, 0.0)
535 for_interior_zones(fn zone
=> update2(cc
,zone
, interior_cc zone
));
536 for_south_zones(fn zone
=> update2(cc
,zone
, reflect_north zone cc
));
537 for_west_zones(fn zone
=> update2(cc
,zone
,reflect_east zone cc
));
538 for_east_zones(fn zone
=> update2(cc
,zone
,reflect_west zone cc
));
539 for_north_zones(fn zone
=> update2(cc
,zone
, reflect_south zone cc
));
543 fun make_sigma(deltat
, rho_prime
, alpha_prime
) =
544 let fun interior_sigma zone
=
545 sub2(rho_prime
,zone
)*sub2(alpha_prime
,zone
)*specific_heat
/ deltat
546 val M
= array2(dimension_interior_zones
, 0.0)
548 (print (Real.toString (sub2(rho_prime
, zone
)) ^
" ");
549 print (Real.toString (sub2(alpha_prime
, zone
)) ^
" ");
550 print (Real.toString specific_heat ^
" ");
551 print (Real.toString deltat ^
"\n");
555 then print ("\t\tmake_sigma:deltat = " ^
Real.toString deltat ^
"\n")
557 (*** for_interior_zones(fn zone
=> update2(M
,zone
, interior_sigma zone
)) **)
558 for_interior_zones(fn zone
=> (update2(M
,zone
, interior_sigma zone
)
559 handle Overflow
=> ohandle zone
));
563 fun make_gamma ((r_prime
,z_prime
), cc
, succeeding
, adjacent
) =
564 let fun interior_gamma zone
=
565 let val r1
= sub2(r_prime
, zone_corner_southeast zone
)
566 val z1
= sub2(z_prime
, zone_corner_southeast zone
)
567 val r2
= sub2(r_prime
, zone_corner_southeast (adjacent zone
))
568 val z2
= sub2(z_prime
, zone_corner_southeast (adjacent zone
))
569 val cross_section
= 0.5*(r1
+r2
)*(pow(r1
- r2
,2)+pow(z1
- z2
,2))
570 val (c1
,c2
) = (sub2(cc
, zone
), sub2(cc
, succeeding zone
))
571 val specific_conductivity
= 2.0 * c1
* c2
/ (c1
+ c2
)
572 in cross_section
* specific_conductivity
574 val M
= array2(dimension_all_zones
, 0.0)
576 for_interior_zones(fn zone
=> update2(M
,zone
,interior_gamma zone
));
580 fun make_ab(theta
, sigma
, Gamma
, preceding
) =
581 let val a
= array2(dimension_all_zones
, 0.0)
582 val b
= array2(dimension_all_zones
, 0.0)
583 fun interior_ab zone
=
584 let val denom
= sub2(sigma
, zone
) + sub2(Gamma
, zone
) +
585 sub2(Gamma
, preceding zone
) *
586 (1.0 - sub2(a
, preceding zone
))
587 val nume1
= sub2(Gamma
,zone
)
588 val nume2
= sub2(Gamma
,preceding zone
)*sub2(b
,preceding zone
) +
589 sub2(sigma
,zone
) * sub2(theta
,zone
)
590 in (nume1
/denom
, nume2
/ denom
)
592 val f
= fn zone
=> update2(b
,zone
,sub2(theta
,zone
))
598 for_interior_zones(fn zone
=> let val ab
= interior_ab zone
599 in update2(a
,zone
,#
1 ab
);
600 update2(b
,zone
,#
2 ab
)
605 fun make_theta (a
, b
, succeeding
, int_zones
) =
606 let val theta
= array2(dimension_all_zones
, constant_heat_source
)
607 fun interior_theta zone
=
608 sub2(a
,zone
) * sub2(theta
,succeeding zone
)+ sub2(b
,zone
)
610 int_zones (fn (k
,l
) => update2(theta
, (k
,l
), interior_theta (k
,l
)));
614 fun compute_heat_conduction(theta_hat
, deltat
, x
', alpha
', rho
') =
615 let val sigma
= make_sigma(deltat
, rho
', alpha
')
616 val _
= if !Control
.trace
then print
"\tdone make_sigma\n" else ()
618 val cc
= make_cc(alpha
', theta_hat
)
619 val _
= if !Control
.trace
then print
"\tdone make_cc\n" else ()
621 val Gamma_k
= make_gamma( x
', cc
, north
, east
)
622 val _
= if !Control
.trace
then print
"\tdone make_gamma\n" else ()
624 val (a_k
,b_k
) = make_ab(theta_hat
, sigma
, Gamma_k
, north
)
625 val _
= if !Control
.trace
then print
"\tdone make_ab\n" else ()
627 val theta_k
= make_theta(a_k
,b_k
,south
,for_north_ward_interior_zones
)
628 val _
= if !Control
.trace
then print
"\tdone make_theta\n" else ()
630 val Gamma_l
= make_gamma(x
', cc
, west
, south
)
631 val _
= if !Control
.trace
then print
"\tdone make_gamma\n" else ()
633 val (a_l
,b_l
) = make_ab(theta_k
, sigma
, Gamma_l
, west
)
634 val _
= if !Control
.trace
then print
"\tdone make_ab\n" else ()
636 val theta_l
= make_theta(a_l
,b_l
,east
,for_west_ward_interior_zones
)
637 val _
= if !Control
.trace
then print
"\tdone make_theta\n" else ()
638 in (theta_l
, Gamma_k
, Gamma_l
)
643 * Final Pressure
and Energy calculation
645 fun make_pressure(rho
', theta
') =
646 let val p
= array2(dimension_all_zones
, 0.0)
647 fun boundary_p(direction
, zone
) =
648 sub1(pbb
, sub2(nbc
, zone
)) +
649 sub1(pb
,sub2(nbc
,zone
)) * sub2(p
, direction zone
)
652 (fn zone
=> update2(p
,zone
,zonal_pressure(sub2(rho
',zone
),
653 sub2(theta
',zone
))));
654 for_south_zones(fn zone
=> update2(p
,zone
,boundary_p(north
,zone
)));
655 for_east_zones(fn zone
=> update2(p
,zone
,boundary_p(west
,zone
)));
656 for_west_zones(fn zone
=> update2(p
,zone
,boundary_p(east
,zone
)));
657 for_north_zones(fn zone
=> update2(p
,zone
,boundary_p(south
,zone
)));
661 fun make_energy(rho
', theta
') =
662 let val epsilon
' = array2(dimension_all_zones
, 0.0)
665 (fn zone
=> update2(epsilon
', zone
, zonal_energy(sub2(rho
',zone
),
666 sub2(theta
',zone
))));
668 (fn zone
=> update2(epsilon
',zone
, reflect_north zone epsilon
'));
670 (fn zone
=> update2(epsilon
',zone
, reflect_east zone epsilon
'));
672 (fn zone
=> update2(epsilon
',zone
, reflect_west zone epsilon
'));
674 (fn zone
=> update2(epsilon
',zone
, reflect_south zone epsilon
'));
680 * Energy Error
Calculation (page
20)
683 fun compute_energy_error ((u
',w
'),(r
',z
'),p
',q
',epsilon
',theta
',rho
',alpha
',
684 Gamma_k
,Gamma_l
,deltat
) =
685 let fun mass zone
= sub2(rho
',zone
) * sub2(alpha
',zone
):real
686 val internal_energy
=
687 sum_list (map_interior_zones (fn z
=> sub2(epsilon
',z
)*(mass z
)))
689 let val average_mass
= 0.25*((mass (zone_A node
)) +
690 (mass (zone_B node
)) +
691 (mass (zone_C node
)) +
692 (mass (zone_D node
)))
693 val v_square
= pow(sub2(u
',node
),2) + pow(sub2(w
',node
),2)
694 in 0.5 * average_mass
* v_square
696 val kinetic_energy
= sum_list (map_interior_nodes kinetic
)
697 fun work_done (node1
, node2
) =
698 let val (r1
, r2
) = (sub2(r
',node1
), sub2(r
',node2
))
699 val (z1
, z2
) = (sub2(z
',node1
), sub2(z
',node2
))
700 val (u1
, u2
) = (sub2(p
',node1
), sub2(p
',node2
))
701 val (w1
, w2
) = (sub2(z
',node1
), sub2(z
',node2
))
702 val (p1
, p2
) = (sub2(p
',node1
), sub2(p
',node2
))
703 val (q1
, q2
) = (sub2(q
',node1
), sub2(q
',node2
))
704 val force
= 0.5*(p1
+p2
+q1
+q2
)
705 val radius
= 0.5* (r1
+r2
)
706 val area
= 0.5* ((r1
-r2
)*(u1
-u2
) - (z1
-z2
)*(w1
-w2
))
707 in force
* radius
* area
* deltat
710 fun from(n
,m
) = if n
> m
then [] else n
::from(n
+1,m
)
712 map (fn l
=> (west(kmin
,l
),(kmin
,l
))) (from(lmin
+1,lmax
))
714 map (fn l
=> (west(kmax
,l
),(kmax
,l
))) (from(lmin
+1,lmax
))
716 map (fn k
=> (south(k
,lmax
),(k
,lmax
))) (from(kmin
+1,kmax
))
718 map (fn k
=> (south(k
,lmin
+1),(k
,lmin
+1))) (from(kmin
+1,kmax
))
720 val w1
= sum_list (map work_done north_line
)
721 val w2
= sum_list (map work_done south_line
)
722 val w3
= sum_list (map work_done east_line
)
723 val w4
= sum_list (map work_done west_line
)
724 val boundary_work
= w1
+ w2
+ w3
+ w4
726 fun heat_flow
Gamma (zone1
,zone2
) =
727 deltat
* sub2(Gamma
, zone1
) * (sub2(theta
',zone1
) - sub2(theta
',zone2
))
731 in map (fn l
=> (north(k
,l
),(k
,l
))) (from(lmin
+1,lmax
))
735 in map (fn l
=> (south(k
,l
),(k
,l
))) (from(lmin
+2,lmax
-1))
739 in map (fn k
=> (east(k
,l
),(k
,l
))) (from(kmin
+2,kmax
))
743 in map (fn k
=> (west(k
,l
),(k
,l
))) (from(kmin
+2,kmax
))
746 val h1
= sum_list (map (heat_flow Gamma_k
) north_flow
)
747 val h2
= sum_list (map (heat_flow Gamma_k
) south_flow
)
748 val h3
= sum_list (map (heat_flow Gamma_l
) east_flow
)
749 val h4
= sum_list (map (heat_flow Gamma_l
) west_flow
)
750 val boundary_heat
= h1
+ h2
+ h3
+ h4
752 internal_energy
+ kinetic_energy
- boundary_heat
- boundary_work
755 fun compute_time_step(d
, theta_hat
, theta
') =
756 let val deltat_courant
=
757 min_list (map_interior_zones (fn zone
=> sub2(d
,zone
)))
759 max_list (map_interior_zones
760 (fn z
=> (abs(sub2(theta_hat
,z
) - sub2(theta
', z
))/
762 val deltat_minimum
= min (deltat_courant
, deltat_conduct
)
763 in min (deltat_maximum
, deltat_minimum
)
767 fun compute_initial_state () =
769 val v
= (all_zero_nodes
, all_zero_nodes
)
770 val x
= let fun interior_position (k
,l
) =
771 let val pi
= 3.1415926535898
772 val rp
= real (lmax
- lmin
)
773 val z1
= real(10 + k
- kmin
)
774 val zz
= (~
0.5 + real(l
- lmin
) / rp
) * pi
775 in (z1
* Math
.cos zz
, z1
* Math
.sin zz
)
777 in make_position_matrix interior_position
780 let val (alpha_prime
,s_prime
) =
781 let val A
= array2(dimension_all_zones
, 0.0)
782 val S
= array2(dimension_all_zones
, 0.0)
783 fun reflect_area_vol f
= (f A
, f S
)
786 let val (a
,s
) = reflect_area_vol(f z
)
792 (fn z
=> let val (a
,s
) = zone_area_vol(x
, z
)
796 for_south_zones (fn z
=> u2 (reflect_north
, z
));
797 for_east_zones (fn z
=> u2 (reflect_west
, z
));
798 for_west_zones (fn z
=> u2 (reflect_east
, z
));
799 for_north_zones (fn z
=> u2 (reflect_south
, z
));
802 in (alpha_prime
,s_prime
)
804 val rho
= let val R
= array2(dimension_all_zones
, 0.0)
805 in for_all_zones (fn z
=> update2(R
,z
,1.4)); R
808 let val T
= array2(dimension_all_zones
, constant_heat_source
)
809 in for_interior_zones(fn z
=> update2(T
,z
,0.0001));
812 val p
= make_pressure(rho
, theta
)
813 val q
= all_zero_zones
814 val epsilon
= make_energy(rho
, theta
)
818 (v
,x
,alpha
,s
,rho
,p
,q
,epsilon
,theta
,deltat
,c
)
822 fun compute_next_state state
=
824 val (v
,x
,alpha
,s
,rho
,p
,q
,epsilon
,theta
,deltat
,c
) = state
825 val v
' = make_velocity (v
, x
, p
, q
, alpha
, rho
, deltat
)
826 val _
= if !Control
.trace
then print
"done make_velocity\n" else ()
828 val x
' = make_position(x
,deltat
,v
')
829 handle Overflow
=>(printarray2 (#
1 v
');
832 val _
= if !Control
.trace
then print
"done make_position\n" else ()
834 val (alpha
',rho
',s
') = make_area_density_volume (rho
, s
, x
')
835 val _
= if !Control
.trace
then print
"done make_area_density_volume\n"
838 val (q
',d
) = make_viscosity (p
, v
', x
', alpha
', rho
')
839 val _
= if !Control
.trace
then print
"done make_viscosity\n" else ()
841 val theta_hat
= make_temperature (p
, epsilon
, rho
, theta
, rho
', q
')
842 val _
= if !Control
.trace
then print
"done make_temperature\n" else ()
844 val (theta
',Gamma_k
,Gamma_l
) =
845 compute_heat_conduction (theta_hat
, deltat
, x
', alpha
', rho
')
846 val _
= if !Control
.trace
then print
"done compute_heat_conduction\n"
849 val p
' = make_pressure(rho
', theta
')
850 val _
= if !Control
.trace
then print
"done make_pressure\n" else ()
852 val epsilon
' = make_energy (rho
', theta
')
853 val _
= if !Control
.trace
then print
"done make_energy\n" else ()
855 val c
' = compute_energy_error (v
', x
', p
', q
', epsilon
', theta
', rho
',
856 alpha
', Gamma_k
, Gamma_l
, deltat
)
857 val _
= if !Control
.trace
then print
"done compute_energy_error\n"
860 val deltat
' = compute_time_step (d
, theta_hat
, theta
')
861 val _
= if !Control
.trace
then print
"done compute_time_step\n\n" else ()
863 (v
',x
',alpha
',s
',rho
',p
',q
', epsilon
',theta
',deltat
',c
')
867 let fun iter (i
,state
) = if i
= 0 then state
869 iter(i
-1, compute_next_state state
))
870 in iter(step_count
, compute_initial_state())
873 fun print_state ((v1
,v2
),(r
,z
),alpha
,s
,rho
,p
,q
,epsilon
,theta
,deltat
,c
) = (
874 print
"Velocity matrices = \n";
875 printarray2 v1
; print
"\n\n";
878 print
"\n\nPosition matrices = \n";
879 printarray2 r
; print
"\n\n";
882 print
"\n\nalpha = \n";
888 print
"\n\nrho = \n";
891 print
"\n\nPressure = \n";
897 print
"\n\nepsilon = \n";
900 print
"\n\ntheta = \n";
903 print ("delatat = " ^
Real.toString (deltat
: real)^
"\n");
904 print ("c = " ^
Real.toString (c
: real) ^
"\n"))
906 fun testit outstrm
= print_state (runit())
909 val (_
, _
, _
, _
, _
, _
, _
, _
, _
, delta
', c
') = runit()
910 val delta
= Real.trunc delta
'
911 val c
= Real.trunc (c
' * 10000.0)
913 if (c
= 6787 andalso delta
= ~
33093)
915 else TextIO.output (TextIO.stdErr
, "*** ERROR ***\n")
930 end; (* functor Simple
*)
931 structure Main
= Simple (val grid_max
=100 val step_count
=1);