HCoop / hcoop / debian / mlton.git / blame
| Commit | Line | Data |
|---|---|---|
| 7f918cf1 CE |
1 | (* From the SML/NJ benchmark suite. *) |
| 2 | fun print _ = () | |
| 3 | structure Control = | |
| 4 | struct | |
| 5 | val trace = ref false | |
| 6 | end; | |
| 7 | structure Array2 : sig | |
| 8 | ||
| 9 | type 'a array2 | |
| 10 | exception Subscript | |
| 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) | |
| 15 | ||
| 16 | end = struct | |
| 17 | ||
| 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 | |
| 22 | else raise Subscript | |
| 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 | |
| 27 | ||
| 28 | end; (* Array2 *) | |
| 29 | signature BMARK = | |
| 30 | sig | |
| 31 | val doit : int -> unit | |
| 32 | val testit : TextIO.outstream -> unit | |
| 33 | end; | |
| 34 | (* Simple | |
| 35 | * error: grid_max < 5 | |
| 36 | *) | |
| 37 | functor Simple(val grid_max: int val step_count: int) : BMARK = | |
| 38 | struct | |
| 39 | ||
| 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)) | |
| 43 | in f2(a,r) | |
| 44 | end) | |
| 45 | ||
| 46 | ||
| 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 | |
| 49 | exception MaxList | |
| 50 | exception MinList | |
| 51 | exception SumList | |
| 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 | |
| 56 | ||
| 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)) | |
| 60 | in f start | |
| 61 | end | |
| 62 | else if endd<=start then | |
| 63 | let fun f x = if x < endd then () else (body x; f(x+delta)) | |
| 64 | in f start | |
| 65 | end | |
| 66 | else () | |
| 67 | fun from(n,m) = if n>m then [] else n::from(n+1,m) | |
| 68 | fun flatten [] = [] | |
| 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) | |
| 76 | fun bounds2(_,b) = b | |
| 77 | fun printarray2 (A as (M:real Array2.array2,((l1,u1),(l2,u2)))) = | |
| 78 | for {from=l1,step=1,to=u1} (fn i => | |
| 79 | (print "["; | |
| 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) | |
| 86 | fun bounds1(_,b) = b | |
| 87 | ||
| 88 | (* | |
| 89 | * Specification of the state variable computation | |
| 90 | *) | |
| 91 | val grid_size = ((2,grid_max), (2,grid_max)) | |
| 92 | ||
| 93 | fun north (k,l) = (k-1,l) | |
| 94 | fun south (k,l) = (k+1,l) | |
| 95 | ||
| 96 | fun east (k,l) = (k,l+1) | |
| 97 | fun west (k,l) = (k,l-1) | |
| 98 | ||
| 99 | val northeast = north o east | |
| 100 | val southeast = south o east | |
| 101 | val northwest = north o west | |
| 102 | val southwest = south o west | |
| 103 | ||
| 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 | |
| 108 | ||
| 109 | fun zone_A(k,l) = (k,l) | |
| 110 | fun zone_B(k,l) = (k+1,l) | |
| 111 | ||
| 112 | fun zone_C(k,l) = (k+1,l+1) | |
| 113 | fun zone_D(k,l) = (k,l+1) | |
| 114 | ||
| 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 | |
| 119 | ||
| 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)) | |
| 125 | ||
| 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)) | |
| 130 | ||
| 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))) | |
| 135 | ||
| 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))) | |
| 140 | ||
| 141 | fun map_interior_nodes f = | |
| 142 | flatten(map (fn k => (map (fn l => f (k,l)) | |
| 143 | (from(lmin,lmax)))) | |
| 144 | (from(kmin,kmax))) | |
| 145 | fun map_interior_zones f = | |
| 146 | flatten(map (fn k => (map (fn l => f (k,l)) | |
| 147 | (from(lmin+1,lmax)))) | |
| 148 | (from(kmin+1,kmax))) | |
| 149 | ||
| 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))) | |
| 156 | ||
| 157 | ||
| 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)) | |
| 162 | ||
| 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 | |
| 168 | ||
| 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)) | |
| 177 | ||
| 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) | |
| 182 | ||
| 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) | |
| 187 | ||
| 188 | ||
| 189 | ||
| 190 | (* | |
| 191 | * Initialization of parameters | |
| 192 | *) | |
| 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 | |
| 198 | end | |
| 199 | val e_coeffs = let val M = array2(((0,2),(0,2)), 0.0) | |
| 200 | in update2(M, (0,1), 0.1); M | |
| 201 | end | |
| 202 | val p_poly = array2(((1,4),(1,5)),p_coeffs) | |
| 203 | ||
| 204 | val e_poly = array2(((1,4),(1,5)), e_coeffs) | |
| 205 | ||
| 206 | val rho_table = let val V = array1((1,3), 0.0) | |
| 207 | in update1(V,2,1.0); | |
| 208 | update1(V,3,100.0); | |
| 209 | V | |
| 210 | end | |
| 211 | val theta_table = let val V = array1((1,4), 0.0) | |
| 212 | in update1(V,2,3.0); | |
| 213 | update1(V,3,300.0); | |
| 214 | update1(V,4,3000.0); | |
| 215 | V | |
| 216 | end | |
| 217 | ||
| 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) | |
| 220 | ||
| 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); | |
| 227 | M | |
| 228 | end | |
| 229 | val pbb = let val A = array1((1,4), 0.0) | |
| 230 | in update1(A,2,6.0); A | |
| 231 | end | |
| 232 | val pb = let val A = array1((1,4), 1.0) | |
| 233 | in update1(A,2,0.0); update1(A,3,0.0); A | |
| 234 | end | |
| 235 | val qb = pb | |
| 236 | ||
| 237 | val all_zero_nodes = array2(dimension_all_nodes, 0.0) | |
| 238 | ||
| 239 | val all_zero_zones = array2(dimension_all_zones, 0.0) | |
| 240 | ||
| 241 | ||
| 242 | (* | |
| 243 | * Positional Coordinates. (page 9-10) | |
| 244 | *) | |
| 245 | ||
| 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 | |
| 255 | in (rb, zb) | |
| 256 | end | |
| 257 | ||
| 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) | |
| 266 | end | |
| 267 | fun u2 (rv,zv) n = (update2(r',n,rv); update2(z',n,zv)) | |
| 268 | in | |
| 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)) | |
| 275 | west_of_north_east; | |
| 276 | u2 (reflect_node(north, northwest, farnorth, west_of_south_east)) | |
| 277 | west_of_south_east; | |
| 278 | u2 (reflect_node(west, northwest, farwest, north_of_south_east)) | |
| 279 | north_of_south_east; | |
| 280 | u2 (reflect_node(east, northeast, fareast, north_of_south_west)) | |
| 281 | north_of_south_west; | |
| 282 | u2 (reflect_node(southwest, west, farwest, north_east_corner)) | |
| 283 | north_east_corner; | |
| 284 | u2 (reflect_node(northwest, west, farwest, south_east_corner)) | |
| 285 | south_east_corner; | |
| 286 | u2 (reflect_node(southeast, south, farsouth, north_west_corner)) | |
| 287 | north_west_corner; | |
| 288 | u2 (reflect_node(northeast, east, fareast, south_west_corner)) | |
| 289 | south_west_corner; | |
| 290 | (r',z') | |
| 291 | end | |
| 292 | ||
| 293 | ||
| 294 | ||
| 295 | (* | |
| 296 | * Physical Properties of a Zone (page 10) | |
| 297 | *) | |
| 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) | |
| 314 | end | |
| 315 | ||
| 316 | (* | |
| 317 | * Velocity (page 8) | |
| 318 | *) | |
| 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)) | |
| 330 | fun velocity 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) | |
| 334 | val u_dot = n1/d | |
| 335 | val w_dot = n2/d | |
| 336 | in (sub2(u,node)+delta_t*u_dot, sub2(w,node)+delta_t*w_dot) | |
| 337 | end | |
| 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); | |
| 342 | update2(W,(k,l),wv) | |
| 343 | end); | |
| 344 | (U,W) | |
| 345 | end | |
| 346 | ||
| 347 | ||
| 348 | ||
| 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 | |
| 354 | end | |
| 355 | ||
| 356 | ||
| 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) | |
| 365 | end | |
| 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); | |
| 369 | update2(s',zone,s); | |
| 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) | |
| 374 | end | |
| 375 | in | |
| 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)); | |
| 381 | (alpha', rho', s') | |
| 382 | end | |
| 383 | ||
| 384 | ||
| 385 | (* | |
| 386 | * Artifical Viscosity (page 11) | |
| 387 | *) | |
| 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))) | |
| 395 | fun lower_del f = | |
| 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 | |
| 408 | val gamma = 1.6 | |
| 409 | val speed_of_sound = gamma*sub2(p,zone)/sub2(rho',zone) | |
| 410 | val ubar = pow(upper_ubar,2) + pow(lower_ubar,2) | |
| 411 | val viscosity = | |
| 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) | |
| 416 | end | |
| 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)) | |
| 423 | in | |
| 424 | for_interior_zones (fn zone => let val (qv,dv) = interior_viscosity zone | |
| 425 | in update2(q',zone,qv); | |
| 426 | update2(d,zone,dv) | |
| 427 | end); | |
| 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)); | |
| 432 | (q', d) | |
| 433 | end | |
| 434 | ||
| 435 | (* | |
| 436 | * Pressure and Energy Polynomial (page 12) | |
| 437 | *) | |
| 438 | ||
| 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) | |
| 444 | in | |
| 445 | if value>sub1(table,high) then high+1 | |
| 446 | else if value <= sub1(table,low) then low | |
| 447 | else search_down high | |
| 448 | end | |
| 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) | |
| 454 | in | |
| 455 | sum_list (map (fn i => sum_list(map (fn j => f (i,j)) (from(0,degree)))) | |
| 456 | (from (0,degree))) | |
| 457 | end | |
| 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) | |
| 462 | end | |
| 463 | ||
| 464 | ||
| 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) | |
| 469 | end | |
| 470 | val dx = 0.000001 | |
| 471 | val tiny = 0.000001 | |
| 472 | ||
| 473 | ||
| 474 | fun newton_raphson (f,x) = | |
| 475 | let fun iter (x,fx) = | |
| 476 | if fx > tiny then | |
| 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) | |
| 481 | end | |
| 482 | else x | |
| 483 | in iter(x, f x) | |
| 484 | end | |
| 485 | ||
| 486 | (* | |
| 487 | * Temperature (page 13-14) | |
| 488 | *) | |
| 489 | ||
| 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 | |
| 510 | in theta_2 | |
| 511 | end | |
| 512 | val M = array2(dimension_all_zones, constant_heat_source) | |
| 513 | in | |
| 514 | for_interior_zones | |
| 515 | (fn zone => update2(M, zone, interior_temperature zone)); | |
| 516 | M | |
| 517 | end | |
| 518 | ||
| 519 | ||
| 520 | (* | |
| 521 | * Heat conduction | |
| 522 | *) | |
| 523 | ||
| 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; | |
| 532 | raise Sqrt) | |
| 533 | val cc = array2(dimension_all_zones, 0.0) | |
| 534 | in | |
| 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)); | |
| 540 | cc | |
| 541 | end | |
| 542 | ||
| 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) | |
| 547 | fun ohandle zone = | |
| 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"); | |
| 552 | raise Overflow) | |
| 553 | ||
| 554 | in if !Control.trace | |
| 555 | then print ("\t\tmake_sigma:deltat = " ^ Real.toString deltat ^ "\n") | |
| 556 | else (); | |
| 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)); | |
| 560 | M | |
| 561 | end | |
| 562 | ||
| 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 | |
| 573 | end | |
| 574 | val M = array2(dimension_all_zones, 0.0) | |
| 575 | in | |
| 576 | for_interior_zones(fn zone => update2(M,zone,interior_gamma zone)); | |
| 577 | M | |
| 578 | end | |
| 579 | ||
| 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) | |
| 591 | end | |
| 592 | val f = fn zone => update2(b,zone,sub2(theta,zone)) | |
| 593 | in | |
| 594 | for_north_zones f; | |
| 595 | for_south_zones f; | |
| 596 | for_west_zones f; | |
| 597 | for_east_zones f; | |
| 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) | |
| 601 | end); | |
| 602 | (a,b) | |
| 603 | end | |
| 604 | ||
| 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) | |
| 609 | in | |
| 610 | int_zones (fn (k,l) => update2(theta, (k,l), interior_theta (k,l))); | |
| 611 | theta | |
| 612 | end | |
| 613 | ||
| 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 () | |
| 617 | ||
| 618 | val cc = make_cc(alpha', theta_hat) | |
| 619 | val _ = if !Control.trace then print "\tdone make_cc\n" else () | |
| 620 | ||
| 621 | val Gamma_k = make_gamma( x', cc, north, east) | |
| 622 | val _ = if !Control.trace then print "\tdone make_gamma\n" else () | |
| 623 | ||
| 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 () | |
| 626 | ||
| 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 () | |
| 629 | ||
| 630 | val Gamma_l = make_gamma(x', cc, west, south) | |
| 631 | val _ = if !Control.trace then print "\tdone make_gamma\n" else () | |
| 632 | ||
| 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 () | |
| 635 | ||
| 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) | |
| 639 | end | |
| 640 | ||
| 641 | ||
| 642 | (* | |
| 643 | * Final Pressure and Energy calculation | |
| 644 | *) | |
| 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) | |
| 650 | in | |
| 651 | for_interior_zones | |
| 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))); | |
| 658 | p | |
| 659 | end | |
| 660 | ||
| 661 | fun make_energy(rho', theta') = | |
| 662 | let val epsilon' = array2(dimension_all_zones, 0.0) | |
| 663 | in | |
| 664 | for_interior_zones | |
| 665 | (fn zone => update2(epsilon', zone, zonal_energy(sub2(rho',zone), | |
| 666 | sub2(theta',zone)))); | |
| 667 | for_south_zones | |
| 668 | (fn zone => update2(epsilon',zone, reflect_north zone epsilon')); | |
| 669 | for_west_zones | |
| 670 | (fn zone => update2(epsilon',zone, reflect_east zone epsilon')); | |
| 671 | for_east_zones | |
| 672 | (fn zone => update2(epsilon',zone, reflect_west zone epsilon')); | |
| 673 | for_north_zones | |
| 674 | (fn zone => update2(epsilon',zone, reflect_south zone epsilon')); | |
| 675 | epsilon' | |
| 676 | end | |
| 677 | ||
| 678 | ||
| 679 | (* | |
| 680 | * Energy Error Calculation (page 20) | |
| 681 | *) | |
| 682 | ||
| 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))) | |
| 688 | fun kinetic node = | |
| 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 | |
| 695 | end | |
| 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 | |
| 708 | end | |
| 709 | ||
| 710 | fun from(n,m) = if n > m then [] else n::from(n+1,m) | |
| 711 | val north_line = | |
| 712 | map (fn l => (west(kmin,l),(kmin,l))) (from(lmin+1,lmax)) | |
| 713 | val south_line = | |
| 714 | map (fn l => (west(kmax,l),(kmax,l))) (from(lmin+1,lmax)) | |
| 715 | val east_line = | |
| 716 | map (fn k => (south(k,lmax),(k,lmax))) (from(kmin+1,kmax)) | |
| 717 | val west_line = | |
| 718 | map (fn k => (south(k,lmin+1),(k,lmin+1))) (from(kmin+1,kmax)) | |
| 719 | ||
| 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 | |
| 725 | ||
| 726 | fun heat_flow Gamma (zone1,zone2) = | |
| 727 | deltat * sub2(Gamma, zone1) * (sub2(theta',zone1) - sub2(theta',zone2)) | |
| 728 | ||
| 729 | val north_flow = | |
| 730 | let val k = kmin+1 | |
| 731 | in map (fn l => (north(k,l),(k,l))) (from(lmin+1,lmax)) | |
| 732 | end | |
| 733 | val south_flow = | |
| 734 | let val k = kmax | |
| 735 | in map (fn l => (south(k,l),(k,l))) (from(lmin+2,lmax-1)) | |
| 736 | end | |
| 737 | val east_flow = | |
| 738 | let val l = lmax | |
| 739 | in map (fn k => (east(k,l),(k,l))) (from(kmin+2,kmax)) | |
| 740 | end | |
| 741 | val west_flow = | |
| 742 | let val l = lmin+1 | |
| 743 | in map (fn k => (west(k,l),(k,l))) (from(kmin+2,kmax)) | |
| 744 | end | |
| 745 | ||
| 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 | |
| 751 | in | |
| 752 | internal_energy + kinetic_energy - boundary_heat - boundary_work | |
| 753 | end | |
| 754 | ||
| 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))) | |
| 758 | val deltat_conduct = | |
| 759 | max_list (map_interior_zones | |
| 760 | (fn z => (abs(sub2(theta_hat,z) - sub2(theta', z))/ | |
| 761 | sub2(theta_hat,z)))) | |
| 762 | val deltat_minimum = min (deltat_courant, deltat_conduct) | |
| 763 | in min (deltat_maximum, deltat_minimum) | |
| 764 | end | |
| 765 | ||
| 766 | ||
| 767 | fun compute_initial_state () = | |
| 768 | let | |
| 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) | |
| 776 | end | |
| 777 | in make_position_matrix interior_position | |
| 778 | end | |
| 779 | val (alpha,s) = | |
| 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) | |
| 784 | ||
| 785 | fun u2 (f,z) = | |
| 786 | let val (a,s) = reflect_area_vol(f z) | |
| 787 | in update2(A,z,a); | |
| 788 | update2(S,z,s) | |
| 789 | end | |
| 790 | in | |
| 791 | for_interior_zones | |
| 792 | (fn z => let val (a,s) = zone_area_vol(x, z) | |
| 793 | in update2(A,z,a); | |
| 794 | update2(S,z,s) | |
| 795 | end); | |
| 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)); | |
| 800 | (A,S) | |
| 801 | end | |
| 802 | in (alpha_prime,s_prime) | |
| 803 | end | |
| 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 | |
| 806 | end | |
| 807 | val theta = | |
| 808 | let val T = array2(dimension_all_zones, constant_heat_source) | |
| 809 | in for_interior_zones(fn z => update2(T,z,0.0001)); | |
| 810 | T | |
| 811 | end | |
| 812 | val p = make_pressure(rho, theta) | |
| 813 | val q = all_zero_zones | |
| 814 | val epsilon = make_energy(rho, theta) | |
| 815 | val deltat = 0.01 | |
| 816 | val c = 0.0 | |
| 817 | in | |
| 818 | (v,x,alpha,s,rho,p,q,epsilon,theta,deltat,c) | |
| 819 | end | |
| 820 | ||
| 821 | ||
| 822 | fun compute_next_state state = | |
| 823 | let | |
| 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 () | |
| 827 | ||
| 828 | val x' = make_position(x,deltat,v') | |
| 829 | handle Overflow =>(printarray2 (#1 v'); | |
| 830 | printarray2 (#2 v'); | |
| 831 | raise Overflow) | |
| 832 | val _ = if !Control.trace then print "done make_position\n" else () | |
| 833 | ||
| 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" | |
| 836 | else () | |
| 837 | ||
| 838 | val (q',d) = make_viscosity (p, v', x', alpha', rho') | |
| 839 | val _ = if !Control.trace then print "done make_viscosity\n" else () | |
| 840 | ||
| 841 | val theta_hat = make_temperature (p, epsilon, rho, theta, rho', q') | |
| 842 | val _ = if !Control.trace then print "done make_temperature\n" else () | |
| 843 | ||
| 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" | |
| 847 | else () | |
| 848 | ||
| 849 | val p' = make_pressure(rho', theta') | |
| 850 | val _ = if !Control.trace then print "done make_pressure\n" else () | |
| 851 | ||
| 852 | val epsilon' = make_energy (rho', theta') | |
| 853 | val _ = if !Control.trace then print "done make_energy\n" else () | |
| 854 | ||
| 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" | |
| 858 | else () | |
| 859 | ||
| 860 | val deltat' = compute_time_step (d, theta_hat, theta') | |
| 861 | val _ = if !Control.trace then print "done compute_time_step\n\n" else () | |
| 862 | in | |
| 863 | (v',x',alpha',s',rho',p',q', epsilon',theta',deltat',c') | |
| 864 | end | |
| 865 | ||
| 866 | fun runit () = | |
| 867 | let fun iter (i,state) = if i = 0 then state | |
| 868 | else (print "."; | |
| 869 | iter(i-1, compute_next_state state)) | |
| 870 | in iter(step_count, compute_initial_state()) | |
| 871 | end | |
| 872 | ||
| 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"; | |
| 876 | printarray2 v2; | |
| 877 | ||
| 878 | print "\n\nPosition matrices = \n"; | |
| 879 | printarray2 r; print "\n\n"; | |
| 880 | printarray2 z; | |
| 881 | ||
| 882 | print "\n\nalpha = \n"; | |
| 883 | printarray2 alpha; | |
| 884 | ||
| 885 | print "\n\ns = \n"; | |
| 886 | printarray2 s; | |
| 887 | ||
| 888 | print "\n\nrho = \n"; | |
| 889 | printarray2 rho; | |
| 890 | ||
| 891 | print "\n\nPressure = \n"; | |
| 892 | printarray2 p; | |
| 893 | ||
| 894 | print "\n\nq = \n"; | |
| 895 | printarray2 q; | |
| 896 | ||
| 897 | print "\n\nepsilon = \n"; | |
| 898 | printarray2 epsilon; | |
| 899 | ||
| 900 | print "\n\ntheta = \n"; | |
| 901 | printarray2 theta; | |
| 902 | ||
| 903 | print ("delatat = " ^ Real.toString (deltat : real)^ "\n"); | |
| 904 | print ("c = " ^ Real.toString (c : real) ^ "\n")) | |
| 905 | ||
| 906 | fun testit outstrm = print_state (runit()) | |
| 907 | ||
| 908 | fun doit () = let | |
| 909 | val (_, _, _, _, _, _, _, _, _, delta', c') = runit() | |
| 910 | val delta = Real.trunc delta' | |
| 911 | val c = Real.trunc (c' * 10000.0) | |
| 912 | in | |
| 913 | if (c = 6787 andalso delta = ~33093) | |
| 914 | then () | |
| 915 | else TextIO.output (TextIO.stdErr, "*** ERROR ***\n") | |
| 916 | end | |
| 917 | ||
| 918 | val doit = | |
| 919 | fn n => | |
| 920 | let | |
| 921 | fun loop n = | |
| 922 | if n = 0 | |
| 923 | then () | |
| 924 | else (doit(); | |
| 925 | loop(n-1)) | |
| 926 | in loop n | |
| 927 | end | |
| 928 | ||
| 929 | ||
| 930 | end; (* functor Simple *) | |
| 931 | structure Main = Simple (val grid_max=100 val step_count=1); |