5 fun digit n
= chr(ord #
"0" + n
)
7 if n
>=0 andalso n
<=9 then digit n
:: acc
8 else digits (n
div 10, digit(n
mod 10) :: acc
)
10 fun int_to_string(n
) =
11 if n
<0 then implode(#
"~"::digits(~n
,[]))
12 else implode(digits(n
,[]))
20 structure Array
= (* Interface
as in SML
/NJ
*)
26 type 'a array
= 'a ref list
33 let fun tab j
= if j
< i
then f j
:: tab (j
+1) else nil
34 in if i
< 0 then raise Size
else (tab
0)
37 fun array (n
, x
) = tabulate
' (n
, fn _
=> ref x
)
39 fun arrayoflist l
= map ref l
41 fun tabulate (n
, f
) = tabulate
' (n
, fn x
=> ref(f x
))
43 fun sub
'(nil
,i
) = raise Subscript
44 | sub
' (a
::r
,i
) = if i
> 0 then sub
' (r
,i
-1)
45 else if i
< 0 then raise Subscript
48 fun op sub (a
, i
) = !(sub
'(a
,i
))
50 fun update (a
, i
, v
) = sub
'(a
, i
) := v
53 |
length (x
::xs
) = 1 + length xs
60 structure List = (* Interface
as in SML
/NJ
*)
64 exception Nth
and NthTail
69 (* this lenght overrides the one
in Array
, but that was also the
70 case in the original lexgen
.sml (see the order
in the
open declaration below
)
74 |
length (x
::xs
) = 1 + length xs
76 fun rev l
= (* linear
-time reversal
of lists
! *)
77 let fun loop([], acc
) = acc
78 |
loop(x
::xs
, acc
) = loop(xs
, x
::acc
)
84 | fold
f (x
::xs
) b
= f(x
,fold f xs b
)
86 fun revfold f
[] b
= b
87 | revfold
f (x
::xs
) b
= revfold f
xs (f(x
,b
))
90 | app
f (x
::xs
) = (f x
; app f xs
)
93 | revapp
f (x
::xs
) = (revapp f xs
; f x
; ())
95 fun nth ([],_
) = raise Nth
97 |
nth (x
::xs
,k
) = nth(xs
,k
-1)
99 fun nthtail ([],_
) = raise NthTail
101 |
nthtail (l
::ls
,n
) = nthtail(ls
,n
-1)
103 fun exists p
[] = false
104 | exists
p (x
::xs
) = p x
orelse exists p xs
121 structure Array2
: sig
124 val array
: (int*int) * '1a
-> '1a array2
125 val sub
: 'a array2
* (int*int) -> 'a
126 val update
: 'a array2
* (int*int) * 'a
-> unit
127 val length
: 'a array2
-> (int*int)
130 type 'a array2
= {size
: (int*int), value
: 'a array
}
131 exception Subscript
= Subscript
132 fun index22 ((i1
:int,i2
:int),(s1
,s2
)) =
133 if i1
>=0 andalso i1
<s1
andalso i2
>=0 andalso i2
<s2
then i1
*s2
+i2
135 fun array22(bnds
as (i1
,i2
), v
) = {size
=bnds
, value
=array(i1
*i2
, v
)}
136 fun op sub22 ({size
,value
}, indx
) = (op sub
) (value
, index22(indx
,size
))
137 fun update22 ({size
=size
,value
=A
},i
,v
) = update(A
,index22(i
,size
),v
)
138 fun length22
{size
=size
,value
=A
} = size
144 * error
: grid_max
< 5
147 functor Simple(val grid_max
: int val step_count
: int) (* : BMARK
*) =
151 (* lars
: unfolding the simple
functor application
*)
161 fun min (x
:real,y
:real) = if x
<y
then x
else y
162 fun max (x
:real,y
:real) = if x
<y
then y
else x
166 fun max_list
[] = raise MaxList | max_list l
= fold max
l (hd l
)
167 fun min_list
[] = raise MinList | min_list l
= fold min
l (hd l
)
168 fun sum_list
[] = raise SumList
169 |
sum_list (l
:real list
) = fold (op +) l
0.0
171 fun for
{from
=start
:int,step
=delta
:int, to
=endd
:int} body
=
172 if delta
>0 andalso endd
>=start
then
173 let fun f x
= if x
> endd
then () else (body x
; f(x
+delta
))
176 else if endd
<=start
then
177 let fun f x
= if x
< endd
then () else (body x
; f(x
+delta
))
181 fun from(n
,m
) = if n
>m
then [] else n
::from(n
+1,m
)
183 |
flatten (x
::xs
) = x @ flatten xs
184 fun pow(x
:real,y
:int) = if y
=0 then 1.0 else x
* pow(x
,y
-1)
185 fun array2(bounds
as ((l1
,u1
),(l2
,u2
)),v
) =
186 (array22((u1
-l1
+1, u2
-l2
+1),v
), bounds
)
187 fun sub2((A
,((lb1
:int,ub1
:int),(lb2
:int,ub2
:int))),(k
,l
)) =
188 sub22(A
, (k
-lb1
, l
-lb2
))
189 fun update2((A
,((lb1
,_
),(lb2
,_
))),(k
,l
), v
) = update22(A
,(k
-lb1
,l
-lb2
),v
)
191 fun printarray2 (A
as (M
:real array2
,((l1
,u1
),(l2
,u2
)))) =
192 for
{from
=l1
,step
=1,to
=u1
} (fn i
=>
194 for
{from
=l2
,step
=1,to
=u2
-1} (fn j
=>
195 print ( (* makestring(sub2(A
,(i
,j
))) ^
*) ", "));
196 print ( (* makestring (sub2(A
,(i
,u2
))) ^
*) "]\n")))
197 fun array1((l
,u
),v
) = (array(u
-l
+1,v
),(l
,u
))
198 fun sub1((A
,(l
:int,u
:int)),i
:int) = (op sub
)(A
,i
-l
)
199 fun update1((A
,(l
,_
)),i
,v
) = update(A
,i
-l
,v
)
203 * Specification
of the state variable computation
205 val grid_size
= ((2,grid_max
), (2,grid_max
))
207 fun north (k
,l
) = (k
-1,l
)
208 fun south (k
,l
) = (k
+1,l
)
210 fun east (k
,l
) = (k
,l
+1)
211 fun west (k
,l
) = (k
,l
-1)
213 val northeast
= north
o east
214 val southeast
= south
o east
215 val northwest
= north
o west
216 val southwest
= south
o west
218 type dir
= int * int -> int * int
220 val farnorth
: dir
= north
o north
221 val farsouth
: dir
= south
o south
222 val fareast
: dir
= east
o east
223 val farwest
: dir
= west
o west
225 fun zone_A(k
,l
) = (k
,l
)
226 fun zone_B(k
,l
) = (k
+1,l
)
228 fun zone_C(k
,l
) = (k
+1,l
+1)
229 fun zone_D(k
,l
) = (k
,l
+1)
231 val zone_corner_northeast
= north
232 val zone_corner_northwest
= northwest
233 fun zone_corner_southeast zone
= zone
234 val zone_corner_southwest
= west
236 val ((kmin
,kmax
),(lmin
,lmax
)) = grid_size
237 val dimension_all_nodes
= ((kmin
-1,kmax
+1),(lmin
-1,lmax
+1))
238 fun for_all_nodes f
=
239 for
{from
=kmin
-1, step
=1, to
=kmax
+1} (fn k
=>
240 for
{from
=lmin
-1, step
=1, to
=lmax
+1} (fn l
=> f k l
))
242 val dimension_interior_nodes
= ((kmin
,kmax
),(lmin
,lmax
))
243 fun for_interior_nodes f
=
244 for
{from
=kmin
, step
=1, to
=kmax
} (fn k
=>
245 for
{from
=lmin
, step
=1, to
=lmax
} (fn l
=> f k l
))
247 val dimension_all_zones
= ((kmin
,kmax
+1),(lmin
,lmax
+1))
248 fun for_all_zones f
=
249 for
{from
=kmin
, step
=1, to
=kmax
+1} (fn k
=>
250 for
{from
=lmin
, step
=1, to
=lmax
+1} (fn l
=> f (k
,l
)))
252 val dimension_interior_zones
= ((kmin
+1,kmax
),(lmin
+1,lmax
))
253 fun for_interior_zones f
=
254 for
{from
=kmin
+1, step
=1, to
=kmax
} (fn k
=>
255 for
{from
=lmin
+1, step
=1, to
=lmax
} (fn l
=> f (k
,l
)))
257 fun map_interior_nodes f
=
258 flatten(map (fn k
=> (map (fn l
=> f (k
,l
))
261 fun map_interior_zones f
=
262 flatten(map (fn k
=> (map (fn l
=> f (k
,l
))
263 (from(lmin
+1,lmax
))))
266 fun for_north_ward_interior_zones f
=
267 for
{from
=kmax
, step
= ~
1, to
=kmin
+1} (fn k
=>
268 for
{from
=lmin
+1, step
=1, to
=lmax
} (fn l
=> f (k
,l
)))
269 fun for_west_ward_interior_zones f
=
270 for
{from
=kmin
+1, step
=1, to
=kmax
} (fn k
=>
271 for
{from
=lmax
, step
= ~
1, to
=lmin
+1} (fn l
=> f (k
,l
)))
274 fun for_north_zones f
= for
{from
=lmin
, step
=1, to
=lmax
+1} (fn l
=> f (kmin
,l
))
275 fun for_south_zones f
= for
{from
=lmin
+1, step
=1, to
=lmax
} (fn l
=> f (kmax
+1,l
))
276 fun for_east_zones f
= for
{from
=kmin
+1, step
=1, to
=kmax
+1}(fn k
=> f (k
,lmax
+1))
277 fun for_west_zones f
= for
{from
=kmin
+1, step
=1, to
=kmax
+1}(fn k
=> f (k
,lmin
))
279 type 'a reflect_dir
= int * int -> {size
: int * int, value
: 'a ref list
}
280 * ((int * int) * (int * int)) -> 'a
282 fun reflect dir node A
= sub2(A
, dir node
)
283 val reflect_north
: real reflect_dir
= reflect north
284 val reflect_south
: real reflect_dir
= reflect south
285 val reflect_east
: real reflect_dir
= reflect east
286 val reflect_west
: real reflect_dir
= reflect west
288 fun for_north_nodes f
=
289 for
{from
=lmin
, step
=1, to
=lmax
-1} (fn l
=> f (kmin
-1,l
))
290 fun for_south_nodes f
=
291 for
{from
=lmin
, step
=1, to
=lmax
-1} (fn l
=> f (kmax
+1,l
))
292 fun for_east_nodes f
=
293 for
{from
=kmin
, step
=1, to
=kmax
-1} (fn k
=> f (k
,lmax
+1))
294 fun for_west_nodes f
=
295 for
{from
=kmin
, step
=1, to
=kmax
-1} (fn k
=> f (k
,lmin
-1))
297 val north_east_corner
= (kmin
-1,lmax
+1)
298 val north_west_corner
= (kmin
-1,lmin
-1)
299 val south_east_corner
= (kmax
+1,lmax
+1)
300 val south_west_corner
= (kmax
+1,lmin
-1)
302 val west_of_north_east
= (kmin
-1, lmax
)
303 val west_of_south_east
= (kmax
+1, lmax
)
304 val north_of_south_east
= (kmax
, lmax
+1)
305 val north_of_south_west
= (kmax
, lmin
-1)
310 * Initialization
of parameters
312 val constant_heat_source
= 0.0
313 val deltat_maximum
= 0.01
314 val specific_heat
= 0.1
315 val p_coeffs
= let val M
= array2(((0,2),(0,2)), 0.0)
316 in update2(M
, (1,1), 0.06698); M
318 val e_coeffs
= let val M
= array2(((0,2),(0,2)), 0.0)
319 in update2(M
, (0,1), 0.1); M
321 val p_poly
= array2(((1,4),(1,5)),p_coeffs
)
323 val e_poly
= array2(((1,4),(1,5)), e_coeffs
)
325 val rho_table
= let val V
= array1((1,3), 0.0)
330 val theta_table
= let val V
= array1((1,4), 0.0)
337 val extract_energy_tables_from_constants
= (e_poly
,2,rho_table
,theta_table
)
338 val extract_pressure_tables_from_constants
= (p_poly
,2,rho_table
,theta_table
)
340 val nbc
= let val M
= array2(dimension_all_zones
, 1)
341 in for
{from
=lmin
+1,step
=1,to
=lmax
} (fn j
=> update2(M
,(kmax
+1, j
),2));
342 update2(M
,(kmin
,lmin
),4);
343 update2(M
,(kmin
,lmax
+1),4);
344 update2(M
,(kmax
+1,lmin
),4);
345 update2(M
,(kmax
+1,lmax
+1),4);
348 val pbb
= let val A
= array1((1,4), 0.0)
349 in update1(A
,2,6.0); A
351 val pb
= let val A
= array1((1,4), 1.0)
352 in update1(A
,2,0.0); update1(A
,3,0.0); A
356 val all_zero_nodes
= array2(dimension_all_nodes
, 0.0)
358 val all_zero_zones
= array2(dimension_all_zones
, 0.0)
362 * Positional Coordinates
. (page
9-10)
365 fun make_position_matrix interior_function
=
366 let val r
' = array2(dimension_all_nodes
, 0.0)
367 val z
' = array2(dimension_all_nodes
, 0.0)
368 fun boundary_position (rx
,zx
,ry
,zy
,ra
,za
) =
369 let val (rax
, zax
) = (ra
- rx
, za
- zx
)
370 val (ryx
, zyx
) = (ry
- rx
, zy
- zx
)
371 val omega
= 2.0*(rax
*ryx
+ zax
*zyx
)/(ryx
*ryx
+ zyx
*zyx
)
372 val rb
= rx
- rax
+ omega
*ryx
373 val zb
= zx
- zax
+ omega
*zyx
377 fun reflect_node (x_dir
, y_dir
, a_dir
, node
) =
378 let val rx
= reflect x_dir node r
'
379 val zx
= reflect x_dir node z
'
380 val ry
= reflect y_dir node r
'
381 val zy
= reflect y_dir node z
'
382 val ra
= reflect a_dir node r
'
383 val za
= reflect a_dir node z
'
384 in boundary_position (rx
, zx
, ry
, zy
, ra
, za
)
386 fun u2 (rv
,zv
) n
= (update2(r
',n
,rv
); update2(z
',n
,zv
))
388 for_interior_nodes (fn k
=> fn l
=> u2 (interior_function (k
,l
)) (k
,l
));
389 for_north_nodes(fn n
=> u2 (reflect_node(south
,southeast
,farsouth
,n
)) n
);
390 for_south_nodes (fn n
=> u2(reflect_node(north
,northeast
,farnorth
,n
)) n
);
391 for_east_nodes (fn n
=> u2(reflect_node(west
, southwest
, farwest
, n
)) n
);
392 for_west_nodes (fn n
=> u2(reflect_node(east
, southeast
, fareast
, n
)) n
);
393 u2 (reflect_node(south
, southwest
, farsouth
, west_of_north_east
))
395 u2 (reflect_node(north
, northwest
, farnorth
, west_of_south_east
))
397 u2 (reflect_node(west
, northwest
, farwest
, north_of_south_east
))
399 u2 (reflect_node(east
, northeast
, fareast
, north_of_south_west
))
401 u2 (reflect_node(southwest
, west
, farwest
, north_east_corner
))
403 u2 (reflect_node(northwest
, west
, farwest
, south_east_corner
))
405 u2 (reflect_node(southeast
, south
, farsouth
, north_west_corner
))
407 u2 (reflect_node(northeast
, east
, fareast
, south_west_corner
))
415 * Physical Properties
of a
Zone (page
10)
417 fun zone_area_vol ((r
,z
), zone
) =
418 let val (r1
,z1
)=(sub2(r
,zone_corner_southwest zone
),
419 sub2(z
,zone_corner_southwest zone
))
420 val (r2
,z2
)=(sub2(r
,zone_corner_southeast zone
),
421 sub2(z
,zone_corner_southeast zone
))
422 val (r3
,z3
)=(sub2(r
,zone_corner_northeast zone
),
423 sub2(z
,zone_corner_northeast zone
))
424 val (r4
,z4
)=(sub2(r
,zone_corner_northwest zone
),
425 sub2(z
,zone_corner_northwest zone
))
426 val area1
= (r2
-r1
)*(z3
-z1
) - (r3
-r2
)*(z3
-z2
)
427 val radius1
= 0.3333 *(r1
+r2
+r3
)
428 val volume1
= area1
* radius1
429 val area2
= (r3
-r1
)*(z4
-z3
) - (r4
-r3
)*(z3
-z1
)
430 val radius2
= 0.3333 *(r1
+r3
+r4
)
431 val volume2
= area2
* radius2
432 in (area1
+area2
, volume1
+volume2
)
438 fun make_velocity((u
,w
),(r
,z
),p
,q
,alpha
,rho
,delta_t
: real) =
439 let fun line_integral (p
,z
,node
) : real =
440 sub2(p
,zone_A node
)*(sub2(z
,west node
) - sub2(z
,north node
)) +
441 sub2(p
,zone_B node
)*(sub2(z
,south node
) - sub2(z
,west node
)) +
442 sub2(p
,zone_C node
)*(sub2(z
,east node
) - sub2(z
,south node
)) +
443 sub2(p
,zone_D node
)*(sub2(z
,north node
) - sub2(z
,east node
))
444 fun regional_mass node
=
445 0.5 * (sub2(rho
, zone_A node
)*sub2(alpha
,zone_A node
) +
446 sub2(rho
, zone_B node
)*sub2(alpha
,zone_B node
) +
447 sub2(rho
, zone_C node
)*sub2(alpha
,zone_C node
) +
448 sub2(rho
, zone_D node
)*sub2(alpha
,zone_D node
))
450 let val d
= regional_mass node
451 val n1
= ~
(line_integral(p
,z
,node
)) - line_integral(q
,z
,node
)
452 val n2
= line_integral(p
,r
,node
) + line_integral(q
,r
,node
)
455 in (sub2(u
,node
)+delta_t
*u_dot
, sub2(w
,node
)+delta_t
*w_dot
)
457 val U
= array2(dimension_interior_nodes
,0.0)
458 val W
= array2(dimension_interior_nodes
,0.0)
459 in for_interior_nodes (fn k
=> fn l
=> let val (uv
,wv
) = velocity (k
,l
)
460 in update2(U
,(k
,l
),uv
);
468 fun make_position ((r
,z
),delta_t
:real,(u
',w
')) =
469 let fun interior_position node
=
470 (sub2(r
,node
) + delta_t
*sub2(u
',node
),
471 sub2(z
,node
) + delta_t
*sub2(w
',node
))
472 in make_position_matrix interior_position
476 fun make_area_density_volume(rho
, s
, x
') =
477 let val alpha
' = array2(dimension_all_zones
, 0.0)
478 val s
' = array2(dimension_all_zones
, 0.0)
479 val rho
' = array2(dimension_all_zones
, 0.0)
480 fun interior_area zone
=
481 let val (area
, vol
) = zone_area_vol (x
', zone
)
482 val density
= sub2(rho
,zone
)*sub2(s
,zone
) / vol
483 in (area
,vol
,density
)
485 fun reflect_area_vol_density reflect_function
=
486 (reflect_function alpha
',reflect_function s
',reflect_function rho
')
487 fun update_asr (zone
,(a
,s
,r
)) = (update2(alpha
',zone
,a
);
489 update2(rho
',zone
,r
))
490 fun r_area_vol_den (reflect_dir
,zone
) =
491 let val asr
= reflect_area_vol_density (reflect_dir zone
)
492 in update_asr(zone
, asr
)
495 for_interior_zones (fn zone
=> update_asr(zone
, interior_area zone
));
496 for_south_zones (fn zone
=> r_area_vol_den(reflect_north
, zone
));
497 for_east_zones (fn zone
=> r_area_vol_den(reflect_west
, zone
));
498 for_west_zones (fn zone
=> r_area_vol_den(reflect_east
, zone
));
499 for_north_zones (fn zone
=> r_area_vol_den(reflect_south
, zone
));
505 * Artifical
Viscosity (page
11)
507 fun make_viscosity(p
,(u
',w
'),(r
',z
'), alpha
',rho
') =
508 let fun interior_viscosity zone
=
509 let fun upper_del f
=
510 0.5 * ((sub2(f
,zone_corner_southeast zone
) -
511 sub2(f
,zone_corner_northeast zone
)) +
512 (sub2(f
,zone_corner_southwest zone
) -
513 sub2(f
,zone_corner_northwest zone
)))
515 0.5 * ((sub2(f
,zone_corner_southeast zone
) -
516 sub2(f
,zone_corner_southwest zone
)) +
517 (sub2(f
,zone_corner_northeast zone
) -
518 sub2(f
,zone_corner_northwest zone
)))
519 val xi
= pow(upper_del r
',2) + pow(upper_del z
',2)
520 val eta
= pow(lower_del r
',2) + pow(lower_del z
',2)
521 val upper_disc
= (upper_del r
')*(lower_del w
') -
522 (upper_del z
')*(lower_del u
')
523 val lower_disc
= (upper_del u
')*(lower_del z
') -
524 (upper_del w
') * (lower_del r
')
525 val upper_ubar
= if upper_disc
<0.0 then upper_disc
/xi
else 0.0
526 val lower_ubar
= if lower_disc
<0.0 then lower_disc
/eta
else 0.0
528 val speed_of_sound
= gamma
*sub2(p
,zone
)/sub2(rho
',zone
)
529 val ubar
= pow(upper_ubar
,2) + pow(lower_ubar
,2)
531 sub2(rho
',zone
)*(1.5*ubar
+ 0.5*speed_of_sound
*(Math
.sqrt ubar
))
532 val length
= Math
.sqrt(pow(upper_del r
',2) + pow(lower_del r
',2))
533 val courant_delta
= 0.5* sub2(alpha
',zone
)/(speed_of_sound
*length
)
534 in (viscosity
, courant_delta
)
536 val q
' = array2(dimension_all_zones
, 0.0)
537 val d
= array2(dimension_all_zones
, 0.0)
538 fun reflect_viscosity_cdelta (direction
, zone
) =
539 sub2(q
',direction zone
) * sub1(qb
, sub2(nbc
,zone
))
540 fun do_zones (dir
,zone
) =
541 update2(q
',zone
,reflect_viscosity_cdelta (dir
,zone
))
543 for_interior_zones (fn zone
=> let val (qv
,dv
) = interior_viscosity zone
544 in update2(q
',zone
,qv
);
547 for_south_zones (fn zone
=> do_zones(north
,zone
));
548 for_east_zones (fn zone
=> do_zones(west
,zone
));
549 for_west_zones (fn zone
=> do_zones(east
,zone
));
550 for_north_zones (fn zone
=> do_zones(south
,zone
));
555 * Pressure
and Energy
Polynomial (page
12)
558 fun polynomial(G
,degree
,rho_table
,theta_table
,rho_value
,theta_value
) =
559 let fun table_search (table
, value
: real) =
560 let val (low
, high
) = bounds1 table
562 if value
> sub1(table
,i
-1) then i
563 else search_down (i
-1)
565 if value
>sub1(table
,high
) then high
+1
566 else if value
<= sub1(table
,low
) then low
567 else search_down high
569 val rho_index
= table_search(rho_table
, rho_value
)
570 val theta_index
= table_search(theta_table
, theta_value
)
571 val A
= sub2(G
, (rho_index
, theta_index
))
572 fun from(n
,m
) = if n
>m
then [] else n
::from(n
+1,m
)
573 fun f(i
,j
) = sub2(A
,(i
,j
))*pow(rho_value
,i
)*pow(theta_value
,j
)
575 sum_list (map (fn i
=> sum_list(map (fn j
=> f (i
,j
)) (from(0,degree
))))
578 fun zonal_pressure (rho_value
:real, theta_value
:real) =
580 val (G
,degree
,rho_table
,theta_table
) =
581 extract_pressure_tables_from_constants
582 in polynomial(G
, degree
, rho_table
, theta_table
, rho_value
, theta_value
)
586 fun zonal_energy (rho_value
, theta_value
) =
587 let val (G
, degree
, rho_table
, theta_table
) =
588 extract_energy_tables_from_constants
589 in polynomial(G
, degree
, rho_table
, theta_table
, rho_value
, theta_value
)
595 fun newton_raphson (f
,x
) =
596 let fun iter (x
,fx
) =
598 let val fxdx
= f(x
+dx
)
599 val denom
= fxdx
- fx
600 in if denom
< tiny
then iter(x
,tiny
)
601 else iter(x
-fx
*dx
/denom
, fxdx
)
608 * Temperature (page
13-14)
611 fun make_temperature(p
,epsilon
,rho
,theta
,rho_prime
,q_prime
) =
612 let fun interior_temperature zone
=
613 let val qkl
= sub2(q_prime
,zone
)
614 val rho_kl
= sub2(rho
,zone
)
615 val rho_prime_kl
= sub2(rho_prime
,zone
)
616 val tau_kl
= (1.0 /rho_prime_kl
- 1.0/rho_kl
)
617 fun energy_equation epsilon_kl theta_kl
=
618 epsilon_kl
- zonal_energy(rho_kl
,theta_kl
)
619 val epsilon_0
= sub2(epsilon
,zone
)
620 fun revised_energy pkl
= epsilon_0
- (pkl
+ qkl
) * tau_kl
621 fun revised_temperature epsilon_kl theta_kl
=
622 newton_raphson ((energy_equation epsilon_kl
), theta_kl
)
623 fun revised_pressure theta_kl
= zonal_pressure(rho_kl
, theta_kl
)
624 val p_0
= sub2(p
,zone
)
625 val theta_0
= sub2(theta
,zone
)
626 val epsilon_1
= revised_energy p_0
627 val theta_1
= revised_temperature epsilon_1 theta_0
628 val p_1
= revised_pressure theta_1
629 val epsilon_2
= revised_energy p_1
630 val theta_2
= revised_temperature epsilon_2 theta_1
633 val M
= array2(dimension_all_zones
, constant_heat_source
)
636 (fn zone
=> update2(M
, zone
, interior_temperature zone
));
645 fun make_cc(alpha_prime
, theta_hat
) =
646 let fun interior_cc zone
=
647 (0.0001 * pow(sub2(theta_hat
,zone
),2) *
648 (Math
.sqrt (abs(sub2(theta_hat
,zone
)))) / sub2(alpha_prime
,zone
))
649 handle Sqrt
=> (print ("<real>" (*Real.makestring (sub2(theta_hat
, zone
))*));
650 print ("\nzone =(" (* ^
makestring (#
1 zone
) *) ^
"," ^
651 (* makestring (#
2 zone
) ^
*) ")\n");
652 printarray2 theta_hat
;
654 val cc
= array2(dimension_all_zones
, 0.0)
656 for_interior_zones(fn zone
=> update2(cc
,zone
, interior_cc zone
));
657 for_south_zones(fn zone
=> update2(cc
,zone
, reflect_north zone cc
));
658 for_west_zones(fn zone
=> update2(cc
,zone
,reflect_east zone cc
));
659 for_east_zones(fn zone
=> update2(cc
,zone
,reflect_west zone cc
));
660 for_north_zones(fn zone
=> update2(cc
,zone
, reflect_south zone cc
));
664 fun make_sigma(deltat
, rho_prime
, alpha_prime
) =
665 let fun interior_sigma zone
=
666 sub2(rho_prime
,zone
)*sub2(alpha_prime
,zone
)*specific_heat
/ deltat
667 val M
= array2(dimension_interior_zones
, 0.0)
669 (print ( (* makestring (sub2(rho_prime
, zone
)) ^
*)" ");
670 print ( (* makestring (sub2(alpha_prime
, zone
)) ^
*)" ");
671 print ( (* makestring specific_heat ^
*) " ");
672 print ( (* makestring deltat ^
*) "\n");
676 then print ("\t\tmake_sigma:deltat = " (* ^ makestring deltat
*) ^
"\n")
678 (*** for_interior_zones(fn zone
=> update2(M
,zone
, interior_sigma zone
)) **)
679 for_interior_zones(fn zone
=> (update2(M
,zone
, interior_sigma zone
)
680 handle _
=> (*old
: Overflow
=> *)
685 fun make_gamma ((r_prime
,z_prime
), cc
, succeeding
, adjacent
) =
686 let fun interior_gamma zone
=
687 let val r1
= sub2(r_prime
, zone_corner_southeast zone
)
688 val z1
= sub2(z_prime
, zone_corner_southeast zone
)
689 val r2
= sub2(r_prime
, zone_corner_southeast (adjacent zone
))
690 val z2
= sub2(z_prime
, zone_corner_southeast (adjacent zone
))
691 val cross_section
= 0.5*(r1
+r2
)*(pow(r1
- r2
,2)+pow(z1
- z2
,2))
692 val (c1
,c2
) = (sub2(cc
, zone
), sub2(cc
, succeeding zone
))
693 val specific_conductivity
= 2.0 * c1
* c2
/ (c1
+ c2
)
694 in cross_section
* specific_conductivity
696 val M
= array2(dimension_all_zones
, 0.0)
698 for_interior_zones(fn zone
=> update2(M
,zone
,interior_gamma zone
));
702 fun make_ab(theta
, sigma
, Gamma
, preceding
) =
703 let val a
= array2(dimension_all_zones
, 0.0)
704 val b
= array2(dimension_all_zones
, 0.0)
705 fun interior_ab zone
=
706 let val denom
= sub2(sigma
, zone
) + sub2(Gamma
, zone
) +
707 sub2(Gamma
, preceding zone
) *
708 (1.0 - sub2(a
, preceding zone
))
709 val nume1
= sub2(Gamma
,zone
)
710 val nume2
= sub2(Gamma
,preceding zone
)*sub2(b
,preceding zone
) +
711 sub2(sigma
,zone
) * sub2(theta
,zone
)
712 in (nume1
/denom
, nume2
/ denom
)
714 val f
= fn zone
=> update2(b
,zone
,sub2(theta
,zone
))
720 for_interior_zones(fn zone
=> let val ab
= interior_ab zone
721 in update2(a
,zone
,#
1 ab
);
722 update2(b
,zone
,#
2 ab
)
727 fun make_theta (a
, b
, succeeding
, int_zones
) =
728 let val theta
= array2(dimension_all_zones
, constant_heat_source
)
729 fun interior_theta zone
=
730 sub2(a
,zone
) * sub2(theta
,succeeding zone
)+ sub2(b
,zone
)
732 int_zones (fn (k
,l
) => update2(theta
, (k
,l
), interior_theta (k
,l
)));
736 fun compute_heat_conduction(theta_hat
, deltat
, x
', alpha
', rho
') =
737 let val sigma
= make_sigma(deltat
, rho
', alpha
')
738 val _
= if !trace
then print
"\tdone make_sigma\n" else ()
740 val cc
= make_cc(alpha
', theta_hat
)
741 val _
= if !trace
then print
"\tdone make_cc\n" else ()
743 val Gamma_k
= make_gamma( x
', cc
, north
, east
)
744 val _
= if !trace
then print
"\tdone make_gamma\n" else ()
746 val (a_k
,b_k
) = make_ab(theta_hat
, sigma
, Gamma_k
, north
)
747 val _
= if !trace
then print
"\tdone make_ab\n" else ()
749 val theta_k
= make_theta(a_k
,b_k
,south
,for_north_ward_interior_zones
)
750 val _
= if !trace
then print
"\tdone make_theta\n" else ()
752 val Gamma_l
= make_gamma(x
', cc
, west
, south
)
753 val _
= if !trace
then print
"\tdone make_gamma\n" else ()
755 val (a_l
,b_l
) = make_ab(theta_k
, sigma
, Gamma_l
, west
)
756 val _
= if !trace
then print
"\tdone make_ab\n" else ()
758 val theta_l
= make_theta(a_l
,b_l
,east
,for_west_ward_interior_zones
)
759 val _
= if !trace
then print
"\tdone make_theta\n" else ()
760 in (theta_l
, Gamma_k
, Gamma_l
)
765 * Final Pressure
and Energy calculation
767 fun make_pressure(rho
', theta
') =
768 let val p
= array2(dimension_all_zones
, 0.0)
769 fun boundary_p(direction
, zone
) =
770 sub1(pbb
, sub2(nbc
, zone
)) +
771 sub1(pb
,sub2(nbc
,zone
)) * sub2(p
, direction zone
)
775 update2(p
,zone
,zonal_pressure(sub2(rho
',zone
),
776 sub2(theta
',zone
))));
777 for_south_zones(fn zone
=> update2(p
,zone
,boundary_p(north
,zone
)));
778 for_east_zones(fn zone
=> update2(p
,zone
,boundary_p(west
,zone
)));
779 for_west_zones(fn zone
=> update2(p
,zone
,boundary_p(east
,zone
)));
780 for_north_zones(fn zone
=> update2(p
,zone
,boundary_p(south
,zone
)));
784 fun make_energy(rho
', theta
') =
785 let val epsilon
' = array2(dimension_all_zones
, 0.0)
788 (fn zone
=> update2(epsilon
', zone
, zonal_energy(sub2(rho
',zone
),
789 sub2(theta
',zone
))));
791 (fn zone
=> update2(epsilon
',zone
, reflect_north zone epsilon
'));
793 (fn zone
=> update2(epsilon
',zone
, reflect_east zone epsilon
'));
795 (fn zone
=> update2(epsilon
',zone
, reflect_west zone epsilon
'));
797 (fn zone
=> update2(epsilon
',zone
, reflect_south zone epsilon
'));
803 * Energy Error
Calculation (page
20)
806 fun compute_energy_error ((u
',w
'),(r
',z
'),p
',q
',epsilon
',theta
',rho
',alpha
',
807 Gamma_k
,Gamma_l
,deltat
) =
808 let fun mass zone
= sub2(rho
',zone
) * sub2(alpha
',zone
):real
809 val internal_energy
=
810 sum_list (map_interior_zones (fn z
=> sub2(epsilon
',z
)*(mass z
)))
812 let val average_mass
= 0.25*((mass (zone_A node
)) +
813 (mass (zone_B node
)) +
814 (mass (zone_C node
)) +
815 (mass (zone_D node
)))
816 val v_square
= pow(sub2(u
',node
),2) + pow(sub2(w
',node
),2)
817 in 0.5 * average_mass
* v_square
819 val kinetic_energy
= sum_list (map_interior_nodes kinetic
)
820 fun work_done (node1
, node2
) =
821 let val (r1
, r2
) = (sub2(r
',node1
), sub2(r
',node2
))
822 val (z1
, z2
) = (sub2(z
',node1
), sub2(z
',node2
))
823 val (u1
, u2
) = (sub2(p
',node1
), sub2(p
',node2
))
824 val (w1
, w2
) = (sub2(z
',node1
), sub2(z
',node2
))
825 val (p1
, p2
) = (sub2(p
',node1
), sub2(p
',node2
))
826 val (q1
, q2
) = (sub2(q
',node1
), sub2(q
',node2
))
827 val force
= 0.5*(p1
+p2
+q1
+q2
)
828 val radius
= 0.5* (r1
+r2
)
829 val area
= 0.5* ((r1
-r2
)*(u1
-u2
) - (z1
-z2
)*(w1
-w2
))
830 in force
* radius
* area
* deltat
833 fun from(n
,m
) = if n
> m
then [] else n
::from(n
+1,m
)
835 map (fn l
=> (west(kmin
,l
),(kmin
,l
))) (from(lmin
+1,lmax
))
837 map (fn l
=> (west(kmax
,l
),(kmax
,l
))) (from(lmin
+1,lmax
))
839 map (fn k
=> (south(k
,lmax
),(k
,lmax
))) (from(kmin
+1,kmax
))
841 map (fn k
=> (south(k
,lmin
+1),(k
,lmin
+1))) (from(kmin
+1,kmax
))
843 val w1
= sum_list (map work_done north_line
)
844 val w2
= sum_list (map work_done south_line
)
845 val w3
= sum_list (map work_done east_line
)
846 val w4
= sum_list (map work_done west_line
)
847 val boundary_work
= w1
+ w2
+ w3
+ w4
849 fun heat_flow
Gamma (zone1
,zone2
) =
850 deltat
* sub2(Gamma
, zone1
) * (sub2(theta
',zone1
) - sub2(theta
',zone2
))
854 in map (fn l
=> (north(k
,l
),(k
,l
))) (from(lmin
+1,lmax
))
858 in map (fn l
=> (south(k
,l
),(k
,l
))) (from(lmin
+2,lmax
-1))
862 in map (fn k
=> (east(k
,l
),(k
,l
))) (from(kmin
+2,kmax
))
866 in map (fn k
=> (west(k
,l
),(k
,l
))) (from(kmin
+2,kmax
))
869 val h1
= sum_list (map (heat_flow Gamma_k
) north_flow
)
870 val h2
= sum_list (map (heat_flow Gamma_k
) south_flow
)
871 val h3
= sum_list (map (heat_flow Gamma_l
) east_flow
)
872 val h4
= sum_list (map (heat_flow Gamma_l
) west_flow
)
873 val boundary_heat
= h1
+ h2
+ h3
+ h4
875 internal_energy
+ kinetic_energy
- boundary_heat
- boundary_work
878 fun compute_time_step(d
, theta_hat
, theta
') =
879 let val deltat_courant
=
880 min_list (map_interior_zones (fn zone
=> sub2(d
,zone
)))
882 max_list (map_interior_zones
883 (fn z
=> (abs(sub2(theta_hat
,z
) - sub2(theta
', z
))/
885 val deltat_minimum
= min (deltat_courant
, deltat_conduct
)
886 in min (deltat_maximum
, deltat_minimum
)
890 fun compute_initial_state () =
892 val v
= (all_zero_nodes
, all_zero_nodes
)
893 val x
= let fun interior_position (k
,l
) =
894 let val pi
= 3.1415926535898
895 val rp
= real (lmax
- lmin
)
896 val z1
= real(10 + k
- kmin
)
897 val zz
= (~
0.5 + real(l
- lmin
) / rp
) * pi
898 in (z1
* Math
.cos zz
, z1
* Math
.sin zz
)
900 in make_position_matrix interior_position
903 let val (alpha_prime
,s_prime
) =
904 let val A
= array2(dimension_all_zones
, 0.0)
905 val S
= array2(dimension_all_zones
, 0.0)
906 fun reflect_area_vol f
= (f A
, f S
)
909 let val (a
,s
) = reflect_area_vol(f z
)
915 (fn z
=> let val (a
,s
) = zone_area_vol(x
, z
)
919 for_south_zones (fn z
=> u2 (reflect_north
, z
));
920 for_east_zones (fn z
=> u2 (reflect_west
, z
));
921 for_west_zones (fn z
=> u2 (reflect_east
, z
));
922 for_north_zones (fn z
=> u2 (reflect_south
, z
));
925 in (alpha_prime
,s_prime
)
927 val rho
= let val R
= array2(dimension_all_zones
, 0.0)
928 in for_all_zones (fn z
=> update2(R
,z
,1.4)); R
931 let val T
= array2(dimension_all_zones
, constant_heat_source
)
932 in for_interior_zones(fn z
=> update2(T
,z
,0.0001));
935 val p
= make_pressure(rho
, theta
)
936 val q
= all_zero_zones
937 val epsilon
= make_energy(rho
, theta
)
941 (v
,x
,alpha
,s
,rho
,p
,q
,epsilon
,theta
,deltat
,c
)
945 fun compute_next_state state
=
947 val (v
,x
,alpha
,s
,rho
,p
,q
,epsilon
,theta
,deltat
,c
) = state
948 val v
' = make_velocity (v
, x
, p
, q
, alpha
, rho
, deltat
)
949 val _
= if !trace
then print
"done make_velocity\n" else ()
951 val x
' = make_position(x
,deltat
,v
')
952 handle _
=> ( (* old
: handle Overflow
=> *)
956 val _
= if !trace
then print
"done make_position\n" else ()
958 val (alpha
',rho
',s
') = make_area_density_volume (rho
, s
, x
')
959 val _
= if !trace
then print
"done make_area_density_volume\n"
962 val (q
',d
) = make_viscosity (p
, v
', x
', alpha
', rho
')
963 val _
= if !trace
then print
"done make_viscosity\n" else ()
965 val theta_hat
= make_temperature (p
, epsilon
, rho
, theta
, rho
', q
')
966 val _
= if !trace
then print
"done make_temperature\n" else ()
968 val (theta
',Gamma_k
,Gamma_l
) =
969 compute_heat_conduction (theta_hat
, deltat
, x
', alpha
', rho
')
970 val _
= if !trace
then print
"done compute_heat_conduction\n"
973 val p
' = make_pressure(rho
', theta
')
974 val _
= if !trace
then print
"done make_pressure\n" else ()
976 val epsilon
' = make_energy (rho
', theta
')
977 val _
= if !trace
then print
"done make_energy\n" else ()
979 val c
' = compute_energy_error (v
', x
', p
', q
', epsilon
', theta
', rho
',
980 alpha
', Gamma_k
, Gamma_l
, deltat
)
981 val _
= if !trace
then print
"done compute_energy_error\n"
984 val deltat
' = compute_time_step (d
, theta_hat
, theta
')
985 val _
= if !trace
then print
"done compute_time_step\n\n" else ()
987 (v
',x
',alpha
',s
',rho
',p
',q
', epsilon
',theta
',deltat
',c
')
991 let fun iter (i
,state
) = if i
= 0 then state
993 iter(i
-1, compute_next_state state
))
994 in iter(step_count
, compute_initial_state())
997 fun print_state ((v1
,v2
),(r
,z
),alpha
,s
,rho
,p
,q
,epsilon
,theta
,deltat
,c
) = (
998 print
"Velocity matrices = \n";
999 printarray2 v1
; print
"\n\n";
1002 print
"\n\nPosition matrices = \n";
1003 printarray2 r
; print
"\n\n";
1006 print
"\n\nalpha = \n";
1012 print
"\n\nrho = \n";
1015 print
"\n\nPressure = \n";
1021 print
"\n\nepsilon = \n";
1022 printarray2 epsilon
;
1024 print
"\n\ntheta = \n";
1027 print ("delatat = " (* ^
Real.makestring deltat
*) ^
"\n");
1028 print ("c = " (* ^
Real.makestring c
*) ^
"\n"))
1030 fun testit outstrm
= print_state (runit())
1033 val (_
, _
, _
, _
, _
, _
, _
, _
, _
, delta
', c
') = runit()
1034 val delta
: int = floor (* truncate
*) delta
'
1035 val c
: int = floor (* truncate
*) (c
' * 10000.0)
1036 val _
= print(int_to_string(c
))
1038 val _
= print(int_to_string(delta
))
1041 if (c
= 3072 andalso delta
= ~
61403) (* for grid_max
= 30 *)
1042 (* (c
= 6787 andalso delta
= ~
33093) *)
1044 else print("*** ERROR ***\n")
1045 (*old
: IO
.output (IO
.std_err
, "*** ERROR ***\n") *)
1049 end; (* functor Simple
*)
1051 structure Main
= Simple(val grid_max
=100 val step_count
=1);