regression/kitsimple.sml

   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