Commit | Line | Data |
---|---|---|
7f918cf1 CE |
1 | (*kitsimple.sml*) |
2 | ||
3 | exception Overflow | |
4 | ||
5 | fun digit n = chr(ord #"0" + n) | |
6 | fun digits(n,acc) = | |
7 | if n >=0 andalso n<=9 then digit n:: acc | |
8 | else digits (n div 10, digit(n mod 10) :: acc) | |
9 | ||
10 | fun int_to_string(n) = | |
11 | if n<0 then implode(#"~"::digits(~n,[])) | |
12 | else implode(digits(n,[])) | |
13 | ||
14 | ||
15 | exception Hd | |
16 | fun hd [] = raise Hd | |
17 | | hd (x::xs) = x | |
18 | ||
19 | (* | |
20 | structure Array = (* Interface as in SML/NJ *) | |
21 | struct | |
22 | *) | |
23 | ||
24 | infix sub | |
25 | ||
26 | type 'a array = 'a ref list | |
27 | ||
28 | exception Size | |
29 | ||
30 | exception Subscript | |
31 | ||
32 | fun tabulate' (i,f) = | |
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) | |
35 | end | |
36 | ||
37 | fun array (n, x) = tabulate' (n, fn _ => ref x) | |
38 | ||
39 | fun arrayoflist l = map ref l | |
40 | ||
41 | fun tabulate (n, f) = tabulate' (n, fn x => ref(f x)) | |
42 | ||
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 | |
46 | else a | |
47 | ||
48 | fun op sub (a, i) = !(sub'(a,i)) | |
49 | ||
50 | fun update (a, i, v) = sub'(a, i) := v | |
51 | ||
52 | fun length [] = 0 | |
53 | | length (x::xs) = 1 + length xs | |
54 | ||
55 | (* | |
56 | end (* Array *) | |
57 | *) | |
58 | ||
59 | (* | |
60 | structure List = (* Interface as in SML/NJ *) | |
61 | struct | |
62 | *) | |
63 | ||
64 | exception Nth and NthTail | |
65 | ||
66 | fun null [] = true | |
67 | | null _ = false | |
68 | ||
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) | |
71 | *) | |
72 | ||
73 | fun length [] = 0 | |
74 | | length (x::xs) = 1 + length xs | |
75 | ||
76 | fun rev l = (* linear-time reversal of lists! *) | |
77 | let fun loop([], acc) = acc | |
78 | | loop(x::xs, acc) = loop(xs, x::acc) | |
79 | in | |
80 | loop(l, []) | |
81 | end | |
82 | ||
83 | fun fold f [] b = b | |
84 | | fold f (x::xs) b = f(x,fold f xs b) | |
85 | ||
86 | fun revfold f [] b = b | |
87 | | revfold f (x::xs) b = revfold f xs (f(x,b)) | |
88 | ||
89 | fun app f [] = () | |
90 | | app f (x::xs) = (f x; app f xs) | |
91 | ||
92 | fun revapp f [] = () | |
93 | | revapp f (x::xs) = (revapp f xs; f x; ()) | |
94 | ||
95 | fun nth ([],_) = raise Nth | |
96 | | nth (x::xs,0) = x | |
97 | | nth (x::xs,k) = nth(xs,k-1) | |
98 | ||
99 | fun nthtail ([],_) = raise NthTail | |
100 | | nthtail (l,0) = l | |
101 | | nthtail (l::ls,n) = nthtail(ls,n-1) | |
102 | ||
103 | fun exists p [] = false | |
104 | | exists p (x::xs) = p x orelse exists p xs | |
105 | ||
106 | (* | |
107 | end; (* List *) | |
108 | *) | |
109 | ||
110 | ||
111 | (* | |
112 | structure Control = | |
113 | struct | |
114 | *) | |
115 | val trace = ref true | |
116 | (* | |
117 | end; | |
118 | *) | |
119 | ||
120 | (* | |
121 | structure Array2 : sig | |
122 | type 'a array2 | |
123 | exception Subscript | |
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) | |
128 | end = struct | |
129 | *) | |
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 | |
134 | else raise Subscript | |
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 | |
139 | (* | |
140 | end; (* Array2 *) | |
141 | *) | |
142 | ||
143 | (* Simple | |
144 | * error: grid_max < 5 | |
145 | *) | |
146 | (* | |
147 | functor Simple(val grid_max: int val step_count: int) (* : BMARK *) = | |
148 | struct | |
149 | *) | |
150 | ||
151 | (* lars: unfolding the simple functor application *) | |
152 | val grid_max = 30 | |
153 | val step_count = 1 | |
154 | ||
155 | (* open Array List | |
156 | *) | |
157 | infix 9 sub | |
158 | ||
159 | ||
160 | ||
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 | |
163 | exception MaxList | |
164 | exception MinList | |
165 | exception SumList | |
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 | |
170 | ||
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)) | |
174 | in f start | |
175 | end | |
176 | else if endd<=start then | |
177 | let fun f x = if x < endd then () else (body x; f(x+delta)) | |
178 | in f start | |
179 | end | |
180 | else () | |
181 | fun from(n,m) = if n>m then [] else n::from(n+1,m) | |
182 | fun flatten [] = [] | |
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) | |
190 | fun bounds2(_,b) = b | |
191 | fun printarray2 (A as (M:real array2,((l1,u1),(l2,u2)))) = | |
192 | for {from=l1,step=1,to=u1} (fn i => | |
193 | (print "["; | |
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) | |
200 | fun bounds1(_,b) = b | |
201 | ||
202 | (* | |
203 | * Specification of the state variable computation | |
204 | *) | |
205 | val grid_size = ((2,grid_max), (2,grid_max)) | |
206 | ||
207 | fun north (k,l) = (k-1,l) | |
208 | fun south (k,l) = (k+1,l) | |
209 | ||
210 | fun east (k,l) = (k,l+1) | |
211 | fun west (k,l) = (k,l-1) | |
212 | ||
213 | val northeast = north o east | |
214 | val southeast = south o east | |
215 | val northwest = north o west | |
216 | val southwest = south o west | |
217 | ||
218 | type dir = int * int -> int * int | |
219 | ||
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 | |
224 | ||
225 | fun zone_A(k,l) = (k,l) | |
226 | fun zone_B(k,l) = (k+1,l) | |
227 | ||
228 | fun zone_C(k,l) = (k+1,l+1) | |
229 | fun zone_D(k,l) = (k,l+1) | |
230 | ||
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 | |
235 | ||
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)) | |
241 | ||
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)) | |
246 | ||
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))) | |
251 | ||
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))) | |
256 | ||
257 | fun map_interior_nodes f = | |
258 | flatten(map (fn k => (map (fn l => f (k,l)) | |
259 | (from(lmin,lmax)))) | |
260 | (from(kmin,kmax))) | |
261 | fun map_interior_zones f = | |
262 | flatten(map (fn k => (map (fn l => f (k,l)) | |
263 | (from(lmin+1,lmax)))) | |
264 | (from(kmin+1,kmax))) | |
265 | ||
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))) | |
272 | ||
273 | ||
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)) | |
278 | ||
279 | type 'a reflect_dir = int * int -> {size: int * int, value: 'a ref list} | |
280 | * ((int * int) * (int * int)) -> 'a | |
281 | ||
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 | |
287 | ||
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)) | |
296 | ||
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) | |
301 | ||
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) | |
306 | ||
307 | ||
308 | ||
309 | (* | |
310 | * Initialization of parameters | |
311 | *) | |
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 | |
317 | end | |
318 | val e_coeffs = let val M = array2(((0,2),(0,2)), 0.0) | |
319 | in update2(M, (0,1), 0.1); M | |
320 | end | |
321 | val p_poly = array2(((1,4),(1,5)),p_coeffs) | |
322 | ||
323 | val e_poly = array2(((1,4),(1,5)), e_coeffs) | |
324 | ||
325 | val rho_table = let val V = array1((1,3), 0.0) | |
326 | in update1(V,2,1.0); | |
327 | update1(V,3,100.0); | |
328 | V | |
329 | end | |
330 | val theta_table = let val V = array1((1,4), 0.0) | |
331 | in update1(V,2,3.0); | |
332 | update1(V,3,300.0); | |
333 | update1(V,4,3000.0); | |
334 | V | |
335 | end | |
336 | ||
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) | |
339 | ||
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); | |
346 | M | |
347 | end | |
348 | val pbb = let val A = array1((1,4), 0.0) | |
349 | in update1(A,2,6.0); A | |
350 | end | |
351 | val pb = let val A = array1((1,4), 1.0) | |
352 | in update1(A,2,0.0); update1(A,3,0.0); A | |
353 | end | |
354 | val qb = pb | |
355 | ||
356 | val all_zero_nodes = array2(dimension_all_nodes, 0.0) | |
357 | ||
358 | val all_zero_zones = array2(dimension_all_zones, 0.0) | |
359 | ||
360 | ||
361 | (* | |
362 | * Positional Coordinates. (page 9-10) | |
363 | *) | |
364 | ||
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 | |
374 | in (rb, zb) | |
375 | end | |
376 | ||
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) | |
385 | end | |
386 | fun u2 (rv,zv) n = (update2(r',n,rv); update2(z',n,zv)) | |
387 | in | |
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)) | |
394 | west_of_north_east; | |
395 | u2 (reflect_node(north, northwest, farnorth, west_of_south_east)) | |
396 | west_of_south_east; | |
397 | u2 (reflect_node(west, northwest, farwest, north_of_south_east)) | |
398 | north_of_south_east; | |
399 | u2 (reflect_node(east, northeast, fareast, north_of_south_west)) | |
400 | north_of_south_west; | |
401 | u2 (reflect_node(southwest, west, farwest, north_east_corner)) | |
402 | north_east_corner; | |
403 | u2 (reflect_node(northwest, west, farwest, south_east_corner)) | |
404 | south_east_corner; | |
405 | u2 (reflect_node(southeast, south, farsouth, north_west_corner)) | |
406 | north_west_corner; | |
407 | u2 (reflect_node(northeast, east, fareast, south_west_corner)) | |
408 | south_west_corner; | |
409 | (r',z') | |
410 | end | |
411 | ||
412 | ||
413 | ||
414 | (* | |
415 | * Physical Properties of a Zone (page 10) | |
416 | *) | |
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) | |
433 | end | |
434 | ||
435 | (* | |
436 | * Velocity (page 8) | |
437 | *) | |
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)) | |
449 | fun velocity 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) | |
453 | val u_dot = n1/d | |
454 | val w_dot = n2/d | |
455 | in (sub2(u,node)+delta_t*u_dot, sub2(w,node)+delta_t*w_dot) | |
456 | end | |
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); | |
461 | update2(W,(k,l),wv) | |
462 | end); | |
463 | (U,W) | |
464 | end | |
465 | ||
466 | ||
467 | ||
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 | |
473 | end | |
474 | ||
475 | ||
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) | |
484 | end | |
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); | |
488 | update2(s',zone,s); | |
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) | |
493 | end | |
494 | in | |
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)); | |
500 | (alpha', rho', s') | |
501 | end | |
502 | ||
503 | ||
504 | (* | |
505 | * Artifical Viscosity (page 11) | |
506 | *) | |
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))) | |
514 | fun lower_del f = | |
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 | |
527 | val gamma = 1.6 | |
528 | val speed_of_sound = gamma*sub2(p,zone)/sub2(rho',zone) | |
529 | val ubar = pow(upper_ubar,2) + pow(lower_ubar,2) | |
530 | val viscosity = | |
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) | |
535 | end | |
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)) | |
542 | in | |
543 | for_interior_zones (fn zone => let val (qv,dv) = interior_viscosity zone | |
544 | in update2(q',zone,qv); | |
545 | update2(d,zone,dv) | |
546 | end); | |
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)); | |
551 | (q', d) | |
552 | end | |
553 | ||
554 | (* | |
555 | * Pressure and Energy Polynomial (page 12) | |
556 | *) | |
557 | ||
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 | |
561 | fun search_down i = | |
562 | if value > sub1(table,i-1) then i | |
563 | else search_down (i-1) | |
564 | in | |
565 | if value>sub1(table,high) then high+1 | |
566 | else if value <= sub1(table,low) then low | |
567 | else search_down high | |
568 | end | |
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) | |
574 | in | |
575 | sum_list (map (fn i => sum_list(map (fn j => f (i,j)) (from(0,degree)))) | |
576 | (from (0,degree))) | |
577 | end | |
578 | fun zonal_pressure (rho_value:real, theta_value:real) = | |
579 | let | |
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) | |
583 | end | |
584 | ||
585 | ||
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) | |
590 | end | |
591 | val dx = 0.000001 | |
592 | val tiny = 0.000001 | |
593 | ||
594 | ||
595 | fun newton_raphson (f,x) = | |
596 | let fun iter (x,fx) = | |
597 | if fx > tiny then | |
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) | |
602 | end | |
603 | else x | |
604 | in iter(x, f x) | |
605 | end | |
606 | ||
607 | (* | |
608 | * Temperature (page 13-14) | |
609 | *) | |
610 | ||
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 | |
631 | in theta_2 | |
632 | end | |
633 | val M = array2(dimension_all_zones, constant_heat_source) | |
634 | in | |
635 | for_interior_zones | |
636 | (fn zone => update2(M, zone, interior_temperature zone)); | |
637 | M | |
638 | end | |
639 | ||
640 | ||
641 | (* | |
642 | * Heat conduction | |
643 | *) | |
644 | ||
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; | |
653 | raise Sqrt) | |
654 | val cc = array2(dimension_all_zones, 0.0) | |
655 | in | |
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)); | |
661 | cc | |
662 | end | |
663 | ||
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) | |
668 | fun ohandle zone = | |
669 | (print ( (* makestring (sub2(rho_prime, zone)) ^ *)" "); | |
670 | print ( (* makestring (sub2(alpha_prime, zone)) ^ *)" "); | |
671 | print ( (* makestring specific_heat ^ *) " "); | |
672 | print ( (* makestring deltat ^ *) "\n"); | |
673 | raise Overflow) | |
674 | ||
675 | in if !trace | |
676 | then print ("\t\tmake_sigma:deltat = " (* ^ makestring deltat *) ^ "\n") | |
677 | else (); | |
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 => *) | |
681 | ohandle zone)); | |
682 | M | |
683 | end | |
684 | ||
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 | |
695 | end | |
696 | val M = array2(dimension_all_zones, 0.0) | |
697 | in | |
698 | for_interior_zones(fn zone => update2(M,zone,interior_gamma zone)); | |
699 | M | |
700 | end | |
701 | ||
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) | |
713 | end | |
714 | val f = fn zone => update2(b,zone,sub2(theta,zone)) | |
715 | in | |
716 | for_north_zones f; | |
717 | for_south_zones f; | |
718 | for_west_zones f; | |
719 | for_east_zones f; | |
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) | |
723 | end); | |
724 | (a,b) | |
725 | end | |
726 | ||
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) | |
731 | in | |
732 | int_zones (fn (k,l) => update2(theta, (k,l), interior_theta (k,l))); | |
733 | theta | |
734 | end | |
735 | ||
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 () | |
739 | ||
740 | val cc = make_cc(alpha', theta_hat) | |
741 | val _ = if !trace then print "\tdone make_cc\n" else () | |
742 | ||
743 | val Gamma_k = make_gamma( x', cc, north, east) | |
744 | val _ = if !trace then print "\tdone make_gamma\n" else () | |
745 | ||
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 () | |
748 | ||
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 () | |
751 | ||
752 | val Gamma_l = make_gamma(x', cc, west, south) | |
753 | val _ = if !trace then print "\tdone make_gamma\n" else () | |
754 | ||
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 () | |
757 | ||
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) | |
761 | end | |
762 | ||
763 | ||
764 | (* | |
765 | * Final Pressure and Energy calculation | |
766 | *) | |
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) | |
772 | in | |
773 | for_interior_zones | |
774 | (fn 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))); | |
781 | p | |
782 | end | |
783 | ||
784 | fun make_energy(rho', theta') = | |
785 | let val epsilon' = array2(dimension_all_zones, 0.0) | |
786 | in | |
787 | for_interior_zones | |
788 | (fn zone => update2(epsilon', zone, zonal_energy(sub2(rho',zone), | |
789 | sub2(theta',zone)))); | |
790 | for_south_zones | |
791 | (fn zone => update2(epsilon',zone, reflect_north zone epsilon')); | |
792 | for_west_zones | |
793 | (fn zone => update2(epsilon',zone, reflect_east zone epsilon')); | |
794 | for_east_zones | |
795 | (fn zone => update2(epsilon',zone, reflect_west zone epsilon')); | |
796 | for_north_zones | |
797 | (fn zone => update2(epsilon',zone, reflect_south zone epsilon')); | |
798 | epsilon' | |
799 | end | |
800 | ||
801 | ||
802 | (* | |
803 | * Energy Error Calculation (page 20) | |
804 | *) | |
805 | ||
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))) | |
811 | fun kinetic node = | |
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 | |
818 | end | |
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 | |
831 | end | |
832 | ||
833 | fun from(n,m) = if n > m then [] else n::from(n+1,m) | |
834 | val north_line = | |
835 | map (fn l => (west(kmin,l),(kmin,l))) (from(lmin+1,lmax)) | |
836 | val south_line = | |
837 | map (fn l => (west(kmax,l),(kmax,l))) (from(lmin+1,lmax)) | |
838 | val east_line = | |
839 | map (fn k => (south(k,lmax),(k,lmax))) (from(kmin+1,kmax)) | |
840 | val west_line = | |
841 | map (fn k => (south(k,lmin+1),(k,lmin+1))) (from(kmin+1,kmax)) | |
842 | ||
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 | |
848 | ||
849 | fun heat_flow Gamma (zone1,zone2) = | |
850 | deltat * sub2(Gamma, zone1) * (sub2(theta',zone1) - sub2(theta',zone2)) | |
851 | ||
852 | val north_flow = | |
853 | let val k = kmin+1 | |
854 | in map (fn l => (north(k,l),(k,l))) (from(lmin+1,lmax)) | |
855 | end | |
856 | val south_flow = | |
857 | let val k = kmax | |
858 | in map (fn l => (south(k,l),(k,l))) (from(lmin+2,lmax-1)) | |
859 | end | |
860 | val east_flow = | |
861 | let val l = lmax | |
862 | in map (fn k => (east(k,l),(k,l))) (from(kmin+2,kmax)) | |
863 | end | |
864 | val west_flow = | |
865 | let val l = lmin+1 | |
866 | in map (fn k => (west(k,l),(k,l))) (from(kmin+2,kmax)) | |
867 | end | |
868 | ||
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 | |
874 | in | |
875 | internal_energy + kinetic_energy - boundary_heat - boundary_work | |
876 | end | |
877 | ||
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))) | |
881 | val deltat_conduct = | |
882 | max_list (map_interior_zones | |
883 | (fn z => (abs(sub2(theta_hat,z) - sub2(theta', z))/ | |
884 | sub2(theta_hat,z)))) | |
885 | val deltat_minimum = min (deltat_courant, deltat_conduct) | |
886 | in min (deltat_maximum, deltat_minimum) | |
887 | end | |
888 | ||
889 | ||
890 | fun compute_initial_state () = | |
891 | let | |
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) | |
899 | end | |
900 | in make_position_matrix interior_position | |
901 | end | |
902 | val (alpha,s) = | |
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) | |
907 | ||
908 | fun u2 (f,z) = | |
909 | let val (a,s) = reflect_area_vol(f z) | |
910 | in update2(A,z,a); | |
911 | update2(S,z,s) | |
912 | end | |
913 | in | |
914 | for_interior_zones | |
915 | (fn z => let val (a,s) = zone_area_vol(x, z) | |
916 | in update2(A,z,a); | |
917 | update2(S,z,s) | |
918 | end); | |
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)); | |
923 | (A,S) | |
924 | end | |
925 | in (alpha_prime,s_prime) | |
926 | end | |
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 | |
929 | end | |
930 | val theta = | |
931 | let val T = array2(dimension_all_zones, constant_heat_source) | |
932 | in for_interior_zones(fn z => update2(T,z,0.0001)); | |
933 | T | |
934 | end | |
935 | val p = make_pressure(rho, theta) | |
936 | val q = all_zero_zones | |
937 | val epsilon = make_energy(rho, theta) | |
938 | val deltat = 0.01 | |
939 | val c = 0.0 | |
940 | in | |
941 | (v,x,alpha,s,rho,p,q,epsilon,theta,deltat,c) | |
942 | end | |
943 | ||
944 | ||
945 | fun compute_next_state state = | |
946 | let | |
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 () | |
950 | ||
951 | val x' = make_position(x,deltat,v') | |
952 | handle _ => ( (* old: handle Overflow => *) | |
953 | printarray2 (#1 v'); | |
954 | printarray2 (#2 v'); | |
955 | raise Overflow) | |
956 | val _ = if !trace then print "done make_position\n" else () | |
957 | ||
958 | val (alpha',rho',s') = make_area_density_volume (rho, s , x') | |
959 | val _ = if !trace then print "done make_area_density_volume\n" | |
960 | else () | |
961 | ||
962 | val (q',d) = make_viscosity (p, v', x', alpha', rho') | |
963 | val _ = if !trace then print "done make_viscosity\n" else () | |
964 | ||
965 | val theta_hat = make_temperature (p, epsilon, rho, theta, rho', q') | |
966 | val _ = if !trace then print "done make_temperature\n" else () | |
967 | ||
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" | |
971 | else () | |
972 | ||
973 | val p' = make_pressure(rho', theta') | |
974 | val _ = if !trace then print "done make_pressure\n" else () | |
975 | ||
976 | val epsilon' = make_energy (rho', theta') | |
977 | val _ = if !trace then print "done make_energy\n" else () | |
978 | ||
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" | |
982 | else () | |
983 | ||
984 | val deltat' = compute_time_step (d, theta_hat, theta') | |
985 | val _ = if !trace then print "done compute_time_step\n\n" else () | |
986 | in | |
987 | (v',x',alpha',s',rho',p',q', epsilon',theta',deltat',c') | |
988 | end | |
989 | ||
990 | fun runit () = | |
991 | let fun iter (i,state) = if i = 0 then state | |
992 | else (print "."; | |
993 | iter(i-1, compute_next_state state)) | |
994 | in iter(step_count, compute_initial_state()) | |
995 | end | |
996 | ||
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"; | |
1000 | printarray2 v2; | |
1001 | ||
1002 | print "\n\nPosition matrices = \n"; | |
1003 | printarray2 r; print "\n\n"; | |
1004 | printarray2 z; | |
1005 | ||
1006 | print "\n\nalpha = \n"; | |
1007 | printarray2 alpha; | |
1008 | ||
1009 | print "\n\ns = \n"; | |
1010 | printarray2 s; | |
1011 | ||
1012 | print "\n\nrho = \n"; | |
1013 | printarray2 rho; | |
1014 | ||
1015 | print "\n\nPressure = \n"; | |
1016 | printarray2 p; | |
1017 | ||
1018 | print "\n\nq = \n"; | |
1019 | printarray2 q; | |
1020 | ||
1021 | print "\n\nepsilon = \n"; | |
1022 | printarray2 epsilon; | |
1023 | ||
1024 | print "\n\ntheta = \n"; | |
1025 | printarray2 theta; | |
1026 | ||
1027 | print ("delatat = " (* ^ Real.makestring deltat *) ^ "\n"); | |
1028 | print ("c = " (* ^ Real.makestring c *) ^ "\n")) | |
1029 | ||
1030 | fun testit outstrm = print_state (runit()) | |
1031 | ||
1032 | fun doit () = let | |
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)) | |
1037 | val _ = print("\n") | |
1038 | val _ = print(int_to_string(delta)) | |
1039 | val _ = print("\n") | |
1040 | in | |
1041 | if (c = 3072 andalso delta = ~61403) (* for grid_max = 30 *) | |
1042 | (* (c = 6787 andalso delta = ~33093) *) | |
1043 | then () | |
1044 | else print("*** ERROR ***\n") | |
1045 | (*old : IO.output (IO.std_err, "*** ERROR ***\n") *) | |
1046 | end | |
1047 | ||
1048 | (* | |
1049 | end; (* functor Simple *) | |
1050 | ||
1051 | structure Main = Simple(val grid_max=100 val step_count=1); | |
1052 | *) | |
1053 | ||
1054 | val _ = doit(); | |
1055 |