1 --simple.hs-- 2 -- The Simple code: 3 -- translated to haskell by an id to ph translator, hand annotated for toplevel 4 -- types, constant types and array strictness. Toplevel comments were readded 5 -- by hand. 6 7 -- ORIGINAL COMMENTS 8 -- 9 -- %%% Id Example Program 10 -- %%% Title: SIMPLE (with higher order functions) 11 -- %%% Original Author: Kattamuri Ekanadham <EKNATH@ibm.com> 12 13 -- %%% Author's Comments: 14 -- % Date: 1 Nov 89 13:33:44 EST 15 -- % From: Kattamuri Ekanadham <EKNATH@ibm.com> 16 -- % Subject: simple89 17 -- % To: jamey%au-bon-pain.lcs.mit.edu@relay.cs.net 18 19 -- % RPaul: 20 -- % I compared the equations in CSG 273 with this code for the following 21 -- % computations: velocity, position, area/volume/density , artificial-viscosity 22 -- % and heat-conduction. 23 -- % - I fixed the area computation in zone_area. (equation 24) 24 -- % - both equation 24 and the code are still missing the constant 25 -- % factor of 2pi in the volume computation 26 -- % - the code for artifical viscocity (equation 27) were missing upper_disc^2 27 -- % and lower_disc^2 28 -- % - the code for interior_cc (equation 37) was missing theta_hat^5/2 29 -- % - the initial condition for heat is now 10.0 (was .0001) 30 -- % 31 -- % -------------- 32 -- % Jamey: 33 -- % I finally re-constructed the scenario that worked some time ago. 34 -- % The following is a program with all defs and built in constant of 35 -- % size = 10. Try this. This should run well. 36 -- % If so, you may try replacing defs by defsubsts and try. 37 -- % -eknath. 38 -- % -------------- 39 -- % eknath : 9/15/89 : the following code is taken directly from the simple paper 40 -- % SIMPLE89-DEFS-OPEN : This has no defsubsts and all constants are in line. 41 -- % The last function simple is invoked with just the number of iterations. 42 -- % Ideally what we want is to wrap this whole thing into the function simple and 43 -- % make grid-max as the second argument for it. 44 45 46 module Main where 47 import Array 48 import Ix 49 infixr 1 =: 50 type Assoc a b = (a,b) 51 (=:) = (,) 52 53 main = getContents >>= \ userInput -> runSimple (lines userInput) 54 55 runSimple :: [String] -> IO () 56 runSimple inputLines = 57 readInt "Run_simple : " inputLines >>= \ num1 -> 58 if num1 >= 0 then reportSimple num1 else reportSimple2 (- num1) 59 60 readInt:: String -> [String] -> IO Int 61 readInt prompt inputLines = 62 putStr prompt >> 63 (case inputLines of 64 (l1:rest) -> case (reads l1) of 65 [(x,"")] -> return x 66 _ -> error "Error" 67 _ -> error "Eof Error") 68 69 reportSimple :: Int -> IO () 70 reportSimple iterations = 71 putStrLn 72 (mix "\n" (let {(u,v,r,z,alpha,s,rho,p,q,epsilon,theta,deltat,err) = simple_total iterations} 73 in [show "RESULT u,v,r,z,alpha,s,rho,p,q,epsilon,theta,deltat err", 74 show "<" ++ show u ++ "," ++ show v ++ "," ++ show r ++ "," ++ show z ++ "," 75 ++ show alpha ++ "," ++ show s ++ "," ++ show rho ++ "," ++ show p ++ "," 76 ++ show q ++ ","++ show epsilon ++ ","++ show theta ++ "," 77 ++ show deltat ++ ","++ show err ++ ">" ]) 78 ++ "\n" ++ "done") 79 80 81 reportSimple2 :: Int -> IO () 82 reportSimple2 iterations = 83 putStrLn 84 (mix "\n" (let {((mat0,mat1),(mat2,mat3),mat4,mat5,mat6,mat7,mat8,mat9,mat10,a,b) = simple iterations} 85 in [show "RESULT ",show "U",show mat0,show "V",show mat1,show "R", 86 show mat2,show "Z", show mat3,show "Alpha", show mat4,show "S", 87 show mat5,show "Rho", show mat6,show "P", 88 show mat7,show "Q", show mat8,show "epsilon", show mat9,"theta", 89 show mat10,show a, show b]) 90 ++ "\n" ++ "done") 91 92 93 94 mix s [] = [] 95 mix s [x] = x 96 mix s (x:xs) = x ++ s ++ mix s xs 97 98 99 100 -- pH prelude functions 101 type Zone = (Int,Int) 102 type DblArray = (Array Zone Double) 103 type DblVector = (Array (Int) Double) 104 type ArrayDim = (Zone,Zone) 105 type VectorDim = Zone 106 type State = ((DblArray, DblArray), (DblArray, DblArray), DblArray, DblArray, DblArray, DblArray, DblArray, DblArray, DblArray, Double, Double) 107 type Totals = (Double, Double, Double, Double, Double, Double, Double, Double, Double, Double, Double, Double, Double) 108 109 110 arrays_2 :: (Ix a) => (a, a) -> [Assoc a (b, c)] -> (Array a b, Array a c) 111 arrays_2 b ivs = ((array b (array_unzip ivs fst)), 112 (array b (array_unzip ivs snd))) 113 114 arrays_3 :: (Ix a) => (a, a) -> [Assoc a (b, c, d)] -> (Array a b, Array a c, Array a d) 115 arrays_3 b ivs = let {first (x,y,z) = x; 116 second (x,y,z) = y; 117 third (x,y,z) = z} 118 in (array b (array_unzip ivs first), 119 array b (array_unzip ivs second), 120 array b (array_unzip ivs third)) 121 122 -- Strict Arrays, force evaluation of the elements. 123 strictArray resultArray = seq (map evalValue (elems resultArray)) resultArray 124 where evalValue v = seq v v 125 126 strictArrays_2 (a1, a2) = (strictArray a1, strictArray a2) 127 strictArrays_3 (a1, a2, a3) = (strictArray a1, strictArray a2, strictArray a3) 128 129 130 131 132 array_unzip :: [Assoc a b] -> (b -> c) -> [Assoc a c] 133 array_unzip ivs sel = [ i =: (sel vs) | (i, vs) <- ivs] 134 135 pHbounds :: ((a, b), (c, d)) -> ((a, c), (b, d)) 136 pHbounds ((x_low,x_high),(y_low,y_high)) = ((x_low,y_low),(x_high,y_high)) 137 --- end pH prelude 138 --- end of driver 139 --- 140 --- 141 --- TRANSLATED CODE 142 --- 143 {- # INCLUDE "list-library"# -} 144 145 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 146 -- %% Simple (toplevel function) 147 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 148 149 simple :: Int -> State 150 simple step_count = let { 151 (totals,state) = simple_loop (compute_initial_state ()) step_count 152 } in state 153 154 simple_total :: Int -> Totals 155 simple_total step_count = 156 let { 157 (totals,state) = simple_loop (compute_initial_state ()) step_count 158 } in totals 159 160 simple_loop :: (Totals,State) -> Int -> (Totals,State) 161 simple_loop result@(totals,state) step = 162 if step < 1 then result else 163 seq (total_total totals) (simple_loop (compute_next_state state)) 164 (step - (1::Int)) 165 166 total_total (v1, v2, x1, x2, alpha, s, rho, 167 p, q, epsilon, theta, deltat, c) = 168 v1 + v2 + x1 + x2 + alpha + s + rho + 169 p + q + epsilon + theta + deltat + c 170 171 total :: (Ix a) => (Array a Double) -> Double 172 total a = foldl (+) (0.0 :: Double) (elems a) 173 174 total_state :: State -> Totals 175 total_state state = 176 let { 177 ((v1, v2), (x1, x2), alpha, s, rho, p, q, epsilon, theta, deltat, c) = 178 state ; 179 v1' = total v1 ; 180 v2' = total v2 ; 181 x1' = total x1 ; 182 x2' = total x2 ; 183 alpha' = total alpha ; 184 s' = total s ; 185 rho' = total rho ; 186 p' = total p ; 187 q' = total q ; 188 epsilon' = total epsilon ; 189 theta' = total theta ; 190 state' = (v1', v2', x1', x2', alpha', s', rho', 191 p', q', epsilon', theta', deltat, c) 192 } in state' 193 194 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 195 -- %% Compute_next_state 196 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 197 198 compute_next_state :: State -> (Totals,State) 199 compute_next_state state = 200 let { 201 state' = (v', x', alpha', s', rho', p', q', 202 epsilon', theta', deltat', c') ; 203 (v, x, alpha, s, rho, p, q, epsilon, theta, deltat, c) = 204 state ; 205 v' = make_velocity v x p q alpha rho deltat ; 206 x' = make_position x deltat v' ; 207 (alpha', s', rho') = make_area_density_volume 208 rho s x' ; 209 (q', d) = make_viscosity p v' x' alpha' rho' ; 210 theta_hat = make_temperature p epsilon rho theta 211 rho' q' ; 212 p' = make_pressure rho' theta' ; 213 epsilon' = make_energy rho' theta' ; 214 (theta', gamma_k, gamma_l) = 215 compute_heat_conduction theta_hat deltat x' alpha' 216 rho' ; 217 c' = compute_energy_error v' x' p' q' epsilon' 218 theta' rho' alpha' gamma_k gamma_l deltat ; 219 (deltat', delta_t_c, delta_t_h) = 220 compute_time_step d theta_hat theta' 221 } in (total_state state', state') 222 223 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 224 -- %% Compute_initial_state 225 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 226 227 compute_initial_state :: a -> (Totals,State) 228 compute_initial_state _ = 229 let { 230 state = (v, x, alpha, s, rho, p, q, epsilon, 231 theta, deltat, c) ; 232 v = (all_zero_nodes (), all_zero_nodes ()) ; 233 x = 234 let { 235 interior_position (k, l) = 236 let { 237 rp = fromIntegral (lmax - lmin) ; 238 z1 = fromIntegral ((10 + k) - kmin) ; 239 zz = ((- (0.5::Double)) + (fromIntegral (l - lmin)) 240 / rp) * (pi::Double) 241 } in ((z1::Double) * (cos (zz::Double)), z1 * (sin zz)) 242 } in make_position_matrix interior_position ; 243 (alpha, s) = 244 let { 245 reflect_area_vol reflect_function = 246 (reflect_function alpha', reflect_function s') ; 247 (alpha', s') = (arrays_2 (pHbounds dimension_all_zones) 248 ([zone=: zone_area_vol x zone 249 | zone <- interior_zones]++ 250 [zone=: reflect_area_vol (reflect_south zone) 251 | zone <- north_zones]++ 252 [zone=: reflect_area_vol (reflect_north zone) 253 | zone <- south_zones]++ 254 [zone=: reflect_area_vol (reflect_west zone) 255 | zone <- east_zones]++ 256 [zone=: reflect_area_vol (reflect_east zone) 257 | zone <- west_zones])) 258 } in strictArrays_2 (alpha', s') ; 259 rho = strictArray (array (pHbounds dimension_all_zones) 260 ([zone=: (1.4::Double) | zone <- all_zones])) ; 261 p = make_pressure rho theta ; 262 q = all_zero_zones () ; 263 epsilon = make_energy rho theta ; 264 theta = strictArray (array (pHbounds dimension_all_zones) 265 ([zone=: (10.0::Double) | zone <- interior_zones]++ 266 [zone=: constant_heat_source 267 | zone <- north_zones]++ 268 [zone=: constant_heat_source 269 | zone <- south_zones]++ 270 [zone=: constant_heat_source 271 | zone <- east_zones]++ 272 [zone=: constant_heat_source 273 | zone <- west_zones])) ; 274 deltat = (0.01::Double) ; 275 c = (0.0 :: Double) 276 } in (total_state state, state) 277 278 line_integral :: DblArray -> DblArray -> (Int,Int) -> Double 279 line_integral p z node = 280 (((p!zone_a node) * ((z!west node) - (z!north node)) 281 + (p!zone_b node) * ((z!south node) 282 - (z!west node))) + (p!zone_c node) 283 * ((z!east node) - (z!south node))) 284 + (p!zone_d node) * ((z!north node) 285 - (z!east node)) 286 287 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 288 -- %% Make_velocity 289 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 290 291 292 make_velocity :: (DblArray, DblArray) -> (DblArray, DblArray) -> DblArray -> DblArray -> (DblArray) -> (DblArray) -> Double -> (DblArray, DblArray) 293 make_velocity (u, w) (r, z) p q alpha rho deltat = 294 295 let { 296 regional_mass node = (0.5::Double) * ((((rho!zone_a node) * (alpha!zone_a node) 297 + (rho!zone_b node) * (alpha!zone_b node)) 298 + (rho!zone_c node) * (alpha!zone_c node)) 299 + (rho!zone_d node) * (alpha!zone_d node)) ; 300 velocity node = 301 let { 302 d = regional_mass node ; 303 n1 = - (line_integral (p :: DblArray) (z :: DblArray) (node :: (Int,Int))) 304 - line_integral q z node ; 305 n2 = line_integral (p :: DblArray) (r :: DblArray) (node :: (Int,Int)) 306 + line_integral (q :: DblArray) r node ; 307 u_dot = (n1 :: Double) / (d :: Double) ; 308 w_dot = (n2 :: Double) / (d :: Double) 309 } in ((u!node) + deltat * u_dot, 310 (w!node) + deltat * w_dot) 311 } in strictArrays_2 (arrays_2 (pHbounds dimension_interior_nodes) 312 ([node=: velocity node | node <- interior_nodes])) 313 314 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 315 -- %% Make_position 316 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 317 318 make_position :: (DblArray, DblArray) -> Double -> (Array (Int,Int) Double, DblArray) -> (DblArray, DblArray) 319 make_position (r, z) deltat (u', w') = 320 let { 321 interior_position node = 322 ((r!node) + deltat * (u'!node), (z!node) + deltat * (w'!node)) 323 } in make_position_matrix interior_position 324 325 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 326 -- %% Make_area_density_volume 327 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 328 329 make_area_density_volume :: (DblArray) -> (DblArray) -> (DblArray, DblArray) -> (DblArray, DblArray, DblArray) 330 make_area_density_volume rho s x' = 331 let { 332 interior_area zone = 333 let { 334 (area, vol) = zone_area_vol x' zone ; 335 density = ((rho!zone) * (s!zone)) / vol 336 } in (area, vol, density) ; 337 reflect_area_vol_density reflect_function = 338 (reflect_function alpha', 339 reflect_function s', reflect_function rho') ; 340 (alpha', s', rho') = (arrays_3 (pHbounds dimension_all_zones) 341 ([zone=: interior_area zone 342 | zone <- interior_zones]++ 343 [zone=: reflect_area_vol_density 344 (reflect_south zone) 345 | zone <- north_zones]++ 346 [zone=: reflect_area_vol_density 347 (reflect_north zone) 348 | zone <- south_zones]++ 349 [zone=: reflect_area_vol_density 350 (reflect_west zone) 351 | zone <- east_zones]++ 352 [zone=: reflect_area_vol_density 353 (reflect_east zone) 354 | zone <- west_zones])) 355 } in strictArrays_3 (alpha', s', rho') 356 357 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 358 -- %% Make_viscosity 359 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 360 361 362 c0 = 1.22474487 363 364 c1 = 0.50 365 366 make_viscosity :: (DblArray) -> (DblArray, DblArray) -> (DblArray, DblArray) -> (DblArray) -> (DblArray) -> (DblArray, DblArray) 367 make_viscosity p (u', w') (r', z') alpha' rho' = 368 let { 369 interior_viscosity zone = 370 let { 371 upper_del f = (0.5::Double) * (((f!zone_corner_southeast zone) 372 - (f!zone_corner_northeast zone)) 373 + (f!zone_corner_southwest zone) 374 - (f!zone_corner_northwest zone)) ; 375 lower_del f = (0.5::Double) * (((f!zone_corner_southeast zone) 376 - (f!zone_corner_southwest zone)) 377 + (f!zone_corner_northeast zone) 378 - (f!zone_corner_northwest zone)) ; 379 xi = (upper_del r') ^ (2::Int) + (upper_del z') 380 ^ (2::Int) ; 381 eta = (lower_del r') ^ (2::Int) + (lower_del z') 382 ^ (2::Int) ; 383 upper_disc = (upper_del r') * (lower_del w') 384 - (upper_del z') * (lower_del u') ; 385 lower_disc = (upper_del u') * (lower_del z') 386 - (upper_del w') * (lower_del r') ; 387 upper_ubar = if upper_disc < (0.0 :: Double) then 388 upper_disc ^ (2::Int) / xi else 389 (0.0 :: Double) ; 390 lower_ubar = if lower_disc < (0.0 :: Double) then 391 lower_disc ^ (2::Int) / eta else 392 (0.0 :: Double) ; 393 gamma = 1.6 ; 394 speed_of_sound = (gamma * (p!zone)) / (rho'!zone) ; 395 ubar = upper_ubar + lower_ubar ; 396 viscosity = (rho'!zone) * ((c0 ^ (2::Int) / 4.0) * ubar 397 + ((c1 / (2.0::Double)) * speed_of_sound) 398 * (sqrt ubar)) ; 399 length = sqrt ((upper_del r') ^ (2::Int) + (lower_del r') 400 ^ (2::Int)) ; 401 courant_delta = ((0.5::Double) * (alpha'!zone)) 402 / speed_of_sound * length 403 } in (viscosity, courant_delta) ; 404 reflect_viscosity_cdelta direction zone = 405 ((q'!direction zone) * (qb!(nbc!zone)), (0.0 :: Double)) ; 406 (q', d) = (arrays_2 (pHbounds dimension_all_zones) 407 ([zone=: interior_viscosity zone 408 | zone <- interior_zones]++ 409 [zone=: reflect_viscosity_cdelta 410 south zone | zone <- north_zones]++ 411 [zone=: reflect_viscosity_cdelta 412 north zone | zone <- south_zones]++ 413 [zone=: reflect_viscosity_cdelta 414 west zone | zone <- east_zones]++ 415 [zone=: reflect_viscosity_cdelta 416 east zone | zone <- west_zones])) 417 } in strictArrays_2 (q', d) 418 419 polynomial :: (Array (Int, Int) (DblArray)) -> Int -> (Array Int Double) -> (Array Int Double) -> Double -> Double -> Double 420 polynomial g degree rho_table theta_table rho_value theta_value = 421 422 let { 423 table_search table value = 424 let { 425 (low, high) = bounds table ; 426 search_down i = if value > (table!(i - 1)) then 427 i else search_down (i - 1) 428 } in if value > (table!high) then 429 high + 1 else if value <= (table!low) then 430 low else search_down high ; 431 rho_index = table_search rho_table rho_value ; 432 theta_index = table_search theta_table theta_value ; 433 a = g!(rho_index, theta_index) 434 } in sum_list [ ((a!(i, j)) * rho_value ^ (i::Int)) 435 * theta_value ^ j | i <- [0..degree],j <- 436 [0..degree] ] 437 438 zonal_pressure :: Double -> Double -> Double 439 zonal_pressure rho_value theta_value = 440 let { 441 (g, degree, rho_table, theta_table) = extract_pressure_tables_from_constants 442 } in polynomial g degree rho_table theta_table 443 rho_value theta_value 444 445 zonal_energy :: Double -> Double -> Double 446 zonal_energy rho_value theta_value = 447 let { 448 (g, degree, rho_table, theta_table) = extract_energy_tables_from_constants 449 } in polynomial g degree rho_table theta_table 450 rho_value theta_value 451 452 dx :: Double 453 dx = 1.0E-6 ::Double 454 455 tiny :: Double 456 tiny = 1.0E-6 ::Double 457 458 newton_raphson :: (Double -> Double) -> Double -> Double 459 newton_raphson f x = newton_raphson_loop f x (f x) 460 461 newton_raphson_loop :: (Double -> Double) -> Double -> Double -> Double 462 newton_raphson_loop f x fx = if fx > tiny then 463 let { 464 fxdx = f (x + dx) ; 465 denom = fxdx - fx ; 466 (next_x, next_fx) = if denom < tiny then (x, 467 tiny) 468 else (x - (fx * dx) / denom, 469 fxdx) 470 } in newton_raphson_loop f next_x 471 next_fx 472 else x 473 474 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 475 -- %% Make_temperature 476 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 477 478 make_temperature :: (DblArray) -> (DblArray) -> (DblArray) -> (DblArray) -> (DblArray) -> (DblArray) -> DblArray 479 make_temperature p epsilon rho theta rho' q' = 480 let { 481 interior_temperature zone = 482 let { 483 qkl = q'!zone ; 484 rho_kl = rho!zone ; 485 rho'_kl = rho'!zone ; 486 tau_kl = 1.0 / rho'_kl - 1.0 / rho_kl ; 487 energy_equation epsilon_kl theta_kl = 488 epsilon_kl - zonal_energy rho_kl theta_kl ; 489 revised_energy pkl = epsilon_0 - (pkl + qkl) 490 * tau_kl ; 491 revised_temperature epsilon_kl theta_kl = 492 newton_raphson (energy_equation epsilon_kl) 493 theta_kl ; 494 revised_pressure theta_kl = 495 zonal_pressure rho_kl theta_kl ; 496 epsilon_0 = epsilon!zone ; 497 p_0 = p!zone ; 498 theta_0 = theta!zone ; 499 epsilon_1 = revised_energy p_0 ; 500 theta_1 = revised_temperature epsilon_1 501 theta_0 ; 502 p_1 = revised_pressure theta_1 ; 503 epsilon_2 = revised_energy p_1 ; 504 theta_2 = revised_temperature epsilon_2 505 theta_1 506 } in theta_2 507 } in strictArray (array (pHbounds dimension_all_zones) 508 ([zone=: interior_temperature zone 509 | zone <- interior_zones]++ 510 [zone=: constant_heat_source 511 | zone <- north_zones]++ 512 [zone=: constant_heat_source 513 | zone <- south_zones]++ 514 [zone=: constant_heat_source 515 | zone <- east_zones]++ 516 [zone=: constant_heat_source 517 | zone <- west_zones])) 518 519 make_cc :: (DblArray) -> (DblArray) -> DblArray 520 make_cc alpha' theta_hat = 521 let { 522 interior_cc zone = ((0.0001::Double) * ((**) (theta_hat!zone) (2.5::Double))) 523 / (alpha'!zone) ; 524 cc = (array (pHbounds dimension_all_zones) 525 ([zone=: interior_cc zone | zone <- interior_zones]++ 526 [zone=: reflect_south zone cc 527 | zone <- north_zones]++ 528 [zone=: reflect_north zone cc 529 | zone <- south_zones]++ 530 [zone=: reflect_west zone cc 531 | zone <- east_zones]++ 532 [zone=: reflect_east zone cc 533 | zone <- west_zones])) 534 } in strictArray cc 535 536 make_sigma :: Double -> (DblArray) -> (DblArray) -> Array (Int,Int) Double 537 make_sigma deltat rho' alpha' = 538 let { 539 interior_sigma zone = (((rho'!zone) * (alpha'!zone)) 540 * specific_heat) / deltat 541 } in strictArray (array (pHbounds dimension_interior_zones) 542 ([zone=: interior_sigma zone 543 | zone <- interior_zones])) 544 545 make_gamma :: (DblArray, DblArray) -> (DblArray)-> ((Int, Int) -> (Int, Int)) -> ((Int, Int) -> (Int, Int)) -> DblArray 546 make_gamma (r', z') cc succeeding adjacent = 547 let { 548 interior_gamma zone = 549 let { 550 r1 = r'!zone_corner_southeast zone ; 551 z1 = z'!zone_corner_southeast zone ; 552 r2 = r'!zone_corner_southeast (adjacent zone) ; 553 z2 = z'!zone_corner_southeast (adjacent zone) ; 554 cross_section = ((0.5::Double) * (r1 + r2)) * ((r1 - r2) ^ (2::Int) 555 + (z1 - z2) ^ (2::Int)) ; 556 (c1, c2) = (cc!zone, cc!succeeding zone) ; 557 specific_conductivity = 558 (((2.0::Double) * c1) * c2) / (c1 + c2) 559 } in cross_section * specific_conductivity 560 } in strictArray (array (pHbounds dimension_all_zones) 561 ([zone=: interior_gamma zone 562 | zone <- interior_zones]++ 563 [zone=: (0.0 :: Double) | zone <- north_zones]++ 564 [zone=: (0.0 :: Double) | zone <- south_zones]++ 565 [zone=: (0.0 :: Double) | zone <- east_zones]++ 566 [zone=: (0.0 :: Double) | zone <- west_zones])) 567 568 make_ab :: (DblArray) -> (DblArray) -> (DblArray) -> ((Int, Int) -> (Int, Int)) 569 -> (DblArray, DblArray) 570 make_ab theta sigma gamma preceding = 571 let { 572 interior_ab zone = 573 let { 574 denom = ((sigma!zone) + (gamma!zone)) 575 + (gamma!preceding zone) * ((1.0 ::Double) 576 - (a!preceding zone)) ; 577 nume1 = gamma!zone ; 578 nume2 = (gamma!preceding zone) 579 * (b!preceding zone) + (sigma!zone) 580 * (theta!zone); 581 result1 = nume1 / denom; 582 result2 = nume2 / denom 583 } in seq (result1 + result2) (result1,result2) ; 584 (a, b) = (arrays_2 (pHbounds dimension_all_zones) 585 ([zone=: interior_ab zone 586 | zone <- interior_zones]++ 587 [zone=: ((0.0 :: Double), theta!zone) 588 | zone <- north_zones]++ 589 [zone=: ((0.0 :: Double), theta!zone) 590 | zone <- south_zones]++ 591 [zone=: ((0.0 :: Double), theta!zone) 592 | zone <- east_zones]++ 593 [zone=: ((0.0 :: Double), theta!zone) 594 | zone <- west_zones])) 595 } in strictArrays_2 (a, b) 596 597 make_theta :: (DblArray) -> (DblArray) -> ((Int, Int) -> (Int, Int)) -> [(Int, Int)] -> DblArray 598 make_theta a b succeeding int_zones = 599 let { 600 interior_theta zone = (a!zone) * (theta!succeeding zone) 601 + (b!zone) ; 602 theta = (array (pHbounds dimension_all_zones) 603 ([zone=: interior_theta zone 604 | zone <- int_zones]++ 605 [zone=: constant_heat_source 606 | zone <- north_zones]++ 607 [zone=: constant_heat_source 608 | zone <- south_zones]++ 609 [zone=: constant_heat_source 610 | zone <- east_zones]++ 611 [zone=: constant_heat_source 612 | zone <- west_zones])) 613 } in strictArray theta 614 615 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 616 -- %% Make_pressure 617 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 618 619 make_pressure :: (DblArray) -> (DblArray) -> DblArray 620 make_pressure rho' theta' = 621 let { 622 boundary_p direction zone = 623 (pbb!(nbc!zone)) + (pb!(nbc!zone)) * (p!direction zone) ; 624 p = (array (pHbounds dimension_all_zones) 625 ([zone=: zonal_pressure (rho'!zone) (theta'!zone) 626 | zone <- interior_zones]++ 627 [zone=: boundary_p south zone 628 | zone <- north_zones]++ 629 [zone=: boundary_p north zone 630 | zone <- south_zones]++ 631 [zone=: boundary_p west zone 632 | zone <- east_zones]++ 633 [zone=: boundary_p east zone 634 | zone <- west_zones])) 635 } in strictArray p 636 637 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 638 -- %% Make_energy 639 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 640 641 make_energy :: (DblArray) -> (DblArray) -> DblArray 642 make_energy rho' theta' = 643 let { 644 epsilon' = (array (pHbounds dimension_all_zones) 645 ([zone=: zonal_energy (rho'!zone) 646 (theta'!zone) 647 | zone <- interior_zones]++ 648 [zone=: reflect_south zone epsilon' 649 | zone <- north_zones]++ 650 [zone=: reflect_north zone epsilon' 651 | zone <- south_zones]++ 652 [zone=: reflect_west zone epsilon' 653 | zone <- east_zones]++ 654 [zone=: reflect_east zone epsilon' 655 | zone <- west_zones])) 656 } in strictArray epsilon' 657 658 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 659 -- %% Compute_heat_conduction 660 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 661 662 compute_heat_conduction :: (DblArray) -> Double -> (DblArray, DblArray) -> (DblArray) -> (DblArray) -> (DblArray, DblArray, DblArray) 663 compute_heat_conduction theta_hat deltat x' alpha' rho' = 664 665 let { 666 sigma = make_sigma deltat rho' alpha' ; 667 cc = make_cc alpha' theta_hat ; 668 gamma_k = make_gamma x' cc north east ; 669 (a_k, b_k) = make_ab theta_hat sigma gamma_k 670 north ; 671 theta_k = make_theta a_k b_k south north_ward_interior_zones ; 672 gamma_l = make_gamma x' cc west south ; 673 (a_l, b_l) = make_ab theta_k sigma gamma_l west ; 674 theta_l = make_theta a_l b_l east west_ward_interior_zones 675 } in (theta_l, gamma_k, gamma_l) 676 677 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 678 -- %% Compute_energy_error 679 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 680 681 compute_energy_error :: (DblArray, DblArray) -> (DblArray, DblArray) -> (DblArray) -> (DblArray) -> (DblArray) -> (DblArray) -> (DblArray) -> (DblArray) -> (DblArray) -> (DblArray) -> Double -> Double 682 compute_energy_error (u', w') (r', z') p' q' epsilon' theta' rho' alpha' gamma_k gamma_l deltat = 683 684 let { 685 mass zone = (rho'!zone) * (alpha'!zone) ; 686 internal_energy = sum_list [ (epsilon'!zone) * (mass zone) 687 | zone <- interior_zones ] ; 688 kinetic node = 689 let { 690 average_mass = (0.25::Double) * (((mass (zone_a node) + mass (zone_b node)) 691 + mass (zone_c node)) 692 + mass (zone_d node)) ; 693 v_square = (u'!node) ^ (2::Int) + (w'!node) 694 ^ (2::Int) 695 } in ((0.5::Double) * average_mass) * v_square ; 696 kinetic_energy = sum_list [ kinetic node | node <- interior_nodes ] ; 697 work_done (node1, node2) = 698 let { 699 (r1, r2) = (r'!node1, r'!node2) ; 700 (z1, z2) = (z'!node1, z'!node2) ; 701 (u1, u2) = (p'!node1, p'!node2) ; 702 (w1, w2) = (z'!node1, z'!node2) ; 703 (p1, p2) = (p'!node1, p'!node2) ; 704 (q1, q2) = (q'!node1, q'!node2) ; 705 force = (0.5::Double) * (((p1 + p2) + q1) + q2) ; 706 radius = (0.5::Double) * (r1 + r2) ; 707 area = (0.5::Double) * ((r1 - r2) * (u1 - u2) 708 - (z1 - z2) * (w1 - w2)) 709 } in ((force * radius) * area) * deltat ; 710 north_line = [ ((i, j), (k, l)) | k <- [kmin],l <- [lmin + 1..lmax], 711 (i, j) <- [west (k, l)] ] ; 712 south_line = [ ((i, j), (k, l)) | k <- [kmax],l <- [lmin + 1..lmax], 713 (i, j) <- [west (k, l)] ] ; 714 east_line = [ ((i, j), (k, l)) | l <- [lmax],k <- [kmin + 1..kmax], 715 (i, j) <- [south (k, l)] ] ; 716 west_line = [ ((i, j), (k, l)) | l <- [lmin + 1],k <- 717 [kmin + 1..kmax],(i, j) <- 718 [south (k, l)] ] ; 719 w1 = sum_list [ work_done segment | segment <- north_line ] ; 720 w2 = sum_list [ work_done segment | segment <- south_line ] ; 721 w3 = sum_list [ work_done segment | segment <- east_line ] ; 722 w4 = sum_list [ work_done segment | segment <- west_line ] ; 723 boundary_work = ((w1 + w2) + w3) + w4 ; 724 heat_flow gamma (zone1, zone2) = 725 (deltat * (gamma!zone1)) * ((theta'!zone1) - (theta'!zone2)) ; 726 north_flow = [ ((i, j), (k, l)) | k <- [kmin + 1],l <- 727 [lmin + 1..lmax],(i, j) <- 728 [north (k, l)] ] ; 729 south_flow = [ ((i, j), (k, l)) | k <- [kmax],l <- [lmin + 2..lmax - 1], 730 (i, j) <- [south (k, l)] ] ; 731 east_flow = [ ((i, j), (k, l)) | l <- [lmax],k <- [kmin + 2..kmax], 732 (i, j) <- [east (k, l)] ] ; 733 west_flow = [ ((i, j), (k, l)) | l <- [lmin + 1],k <- 734 [kmin + 2..kmax],(i, j) <- 735 [west (k, l)] ] ; 736 h1 = sum_list [ heat_flow gamma_k segment 737 | segment <- north_flow ] ; 738 h2 = sum_list [ heat_flow gamma_k segment 739 | segment <- south_flow ] ; 740 h3 = sum_list [ heat_flow gamma_l segment 741 | segment <- east_flow ] ; 742 h4 = sum_list [ heat_flow gamma_l segment 743 | segment <- west_flow ] ; 744 boundary_heat = ((h1 + h2) + h3) + h4 745 } in ((internal_energy + kinetic_energy) - boundary_heat) 746 - boundary_work 747 748 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 749 -- %% Compute_time_step 750 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 751 752 compute_time_step :: (DblArray) -> (DblArray) -> (DblArray) -> (Double, Double, Double) 753 compute_time_step d theta_hat theta' = 754 let { 755 deltat_courant = min_list [ d!zone | zone <- interior_zones ] ; 756 deltat_conduct = max_list [ (abs ((theta_hat!zone) - (theta'!zone))) 757 / (theta_hat!zone) | zone <- interior_zones ] ; 758 deltat_minimum = min deltat_courant deltat_conduct 759 } in (min deltat_maximum deltat_minimum, deltat_courant, 760 deltat_conduct) 761 762 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 763 -- %% Geometry 764 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 765 766 north :: Zone -> Zone 767 north (k, l) = (k - 1, l) 768 769 south :: Zone -> Zone 770 south (k, l) = (k + 1, l) 771 772 east :: Zone -> Zone 773 east (k, l) = (k, l + 1) 774 775 west :: Zone -> Zone 776 west (k, l) = (k, l - 1) 777 778 northeast :: Zone -> Zone 779 northeast node = north (east node) 780 781 southeast :: Zone -> Zone 782 southeast node = south (east node) 783 784 northwest :: Zone -> Zone 785 northwest node = north (west node) 786 787 southwest :: Zone -> Zone 788 southwest node = south (west node) 789 790 farnorth :: Zone -> Zone 791 farnorth node = north (north node) 792 793 farsouth :: Zone -> Zone 794 farsouth node = south (south node) 795 796 fareast :: Zone -> Zone 797 fareast node = east (east node) 798 799 farwest :: Zone -> Zone 800 farwest node = west (west node) 801 802 zone_a :: Zone -> Zone 803 zone_a (k, l) = (k, l) 804 805 zone_b :: Zone -> Zone 806 zone_b (k, l) = (k + 1, l) 807 808 zone_c :: Zone -> Zone 809 zone_c (k, l) = (k + 1, l + 1) 810 811 zone_d :: Zone -> Zone 812 zone_d (k, l) = (k, l + 1) 813 814 zone_corner_northeast :: Zone -> Zone 815 zone_corner_northeast zone = north zone 816 817 zone_corner_northwest :: Zone -> Zone 818 zone_corner_northwest zone = northwest zone 819 820 zone_corner_southeast :: Zone -> Zone 821 zone_corner_southeast zone = zone 822 823 zone_corner_southwest :: Zone -> Zone 824 zone_corner_southwest zone = west zone 825 826 ((kmin, kmax), (lmin, lmax)) = grid_size :: ArrayDim 827 828 dimension_all_nodes = ((kmin - 1, kmax + 1), 829 (lmin - 1, lmax + 1)) :: ArrayDim 830 831 all_nodes :: [(Int, Int)] 832 all_nodes = [ (k, l) | k <- [kmin - 1..kmax + 1],l <- 833 [lmin - 1..lmax + 1] ] 834 835 dimension_interior_nodes :: ArrayDim 836 dimension_interior_nodes = ((kmin, kmax), (lmin, lmax)) 837 838 interior_nodes :: [(Int, Int)] 839 interior_nodes = [ (k, l) | k <- [kmin..kmax],l <- 840 [lmin..lmax] ] 841 842 dimension_all_zones :: ArrayDim 843 dimension_all_zones = ((kmin, kmax + 1), (lmin, lmax + 1)) 844 845 all_zones :: [(Int, Int)] 846 all_zones = [ (k, l) | k <- [kmin..kmax + 1],l <- 847 [lmin..lmax + 1] ] 848 849 dimension_interior_zones :: ArrayDim 850 dimension_interior_zones = 851 ((kmin + 1, kmax), (lmin + 1, lmax)) 852 853 interior_zones :: [(Int, Int)] 854 interior_zones = [ (k, l) | k <- [kmin + 1..kmax],l <- 855 [lmin + 1..lmax] ] 856 857 north_ward_interior_zones :: [(Int, Int)] 858 north_ward_interior_zones = 859 [ (k, l) | k <- [kmax,(kmax-1)..kmin + 1],l <- [lmin + 1..lmax] ] 860 861 west_ward_interior_zones :: [(Int, Int)] 862 west_ward_interior_zones = 863 [ (k, l) | k <- [kmin + 1..kmax],l <- [lmax,(lmax-1)..lmin + 1] ] 864 865 north_zones :: [(Int, Int)] 866 north_zones = [ (k, l) | k <- [kmin],l <- [lmin..lmax + 1] ] 867 868 south_zones :: [(Int, Int)] 869 south_zones = [ (k, l) | k <- [kmax + 1],l <- 870 [lmin + 1..lmax] ] 871 872 east_zones :: [(Int, Int)] 873 east_zones = [ (k, l) | k <- [kmin + 1..kmax + 1],l <- 874 [lmax + 1] ] 875 876 west_zones :: [(Int, Int)] 877 west_zones = [ (k, l) | k <- [kmin + 1..kmax + 1],l <- 878 [lmin] ] 879 880 reflect :: ((Int, Int) -> (Int, Int)) -> (Int, Int) -> DblArray -> Double 881 reflect dir node a = a!dir node 882 883 reflect_north :: (Int, Int) -> (Array (Int, Int) Double) -> Double 884 reflect_north = reflect north 885 886 reflect_south :: (Int, Int) -> (Array (Int, Int) Double) -> Double 887 reflect_south = reflect south 888 889 reflect_east :: (Int, Int) -> (Array (Int, Int) Double) -> Double 890 reflect_east = reflect east 891 892 reflect_west :: (Int, Int) -> (Array (Int, Int) Double) -> Double 893 reflect_west = reflect west 894 895 north_nodes :: [(Int, Int)] 896 north_nodes = [ (k, l) | k <- [kmin - 1],l <- 897 [lmin..lmax - 1] ] 898 899 south_nodes :: [(Int, Int)] 900 south_nodes = [ (k, l) | k <- [kmax + 1],l <- 901 [lmin..lmax - 1] ] 902 903 east_nodes :: [(Int, Int)] 904 east_nodes = [ (k, l) | k <- [kmin..kmax - 1],l <- 905 [lmax + 1] ] 906 907 west_nodes :: [(Int, Int)] 908 west_nodes = [ (k, l) | k <- [kmin..kmax - 1],l <- 909 [lmin - 1] ] 910 911 north_east_corner :: [(Int, Int)] 912 north_east_corner = [ (k, l) | k <- [kmin - 1],l <- 913 [lmax + 1] ] 914 915 north_west_corner :: [(Int, Int)] 916 north_west_corner = [ (k, l) | k <- [kmin - 1],l <- 917 [lmin - 1] ] 918 919 south_east_corner :: [(Int, Int)] 920 south_east_corner = [ (k, l) | k <- [kmax + 1],l <- 921 [lmax + 1] ] 922 923 south_west_corner :: [(Int, Int)] 924 south_west_corner = [ (k, l) | k <- [kmax + 1],l <- 925 [lmin - 1] ] 926 927 west_of_north_east :: [(Int, Int)] 928 west_of_north_east = [ (k, l) | k <- [kmin - 1],l <- 929 [lmax] ] 930 931 west_of_south_east :: [(Int, Int)] 932 west_of_south_east = [ (k, l) | k <- [kmax + 1],l <- 933 [lmax] ] 934 935 north_of_south_east :: [(Int, Int)] 936 north_of_south_east = [ (k, l) | k <- [kmax],l <- [lmax + 1] ] 937 938 north_of_south_west :: [(Int, Int)] 939 north_of_south_west = [ (k, l) | k <- [kmax],l <- [lmin - 1] ] 940 941 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 942 -- %% Constants 943 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 944 945 -- % eknath : initializations and other additions follow here 946 947 grid_max = 100 :: Int 948 949 grid_size = ((2, grid_max), (2, grid_max)) ::ArrayDim 950 951 constant_heat_source = (0.0 :: Double) 952 953 deltat_maximum = (0.01 :: Double) 954 955 specific_heat = 0.1 ::Double 956 957 p_coeffs :: DblArray 958 p_coeffs = strictArray (array ((0, 0),(2, 2)) 959 ([(1, 1) =: 0.06698]++ 960 [(i, j)=: 0.0 | i <- [0..2],j <- [0..2], not (i == 1 && j == 1)])) 961 962 e_coeffs:: DblArray 963 e_coeffs = strictArray (array ((0, 0),(2, 2)) 964 ([(0, 1) =: 0.1]++ 965 [(i, j)=: 0.0 | i <- [0..2],j <- [0..2], not (i == 0 && j == 1)])) 966 967 p_poly = strictArray (array ((1, 1),(4, 5)) 968 ([(i, j)=: p_coeffs | i <- [1..4],j <- [1..5]])) 969 970 e_poly = strictArray (array ((1, 1),(4, 5)) 971 ([(i, j)=: e_coeffs | i <- [1..4],j <- [1..5]])) 972 973 rho_table:: DblVector 974 rho_table = strictArray (array (1, 3) 975 ([1 =: 0.0]++ 976 [2 =: 1.0]++ 977 [3 =: 100.0])) 978 979 theta_table:: DblVector 980 theta_table = strictArray (array (1, 4) 981 ([1 =: 0.0]++ 982 [2 =: 3.0]++ 983 [3 =: 300.0]++ 984 [4 =: 3000.0])) 985 986 extract_energy_tables_from_constants = 987 (e_poly, 2, rho_table, theta_table) 988 989 extract_pressure_tables_from_constants = 990 (p_poly, 2, rho_table, theta_table) 991 992 nbc :: Array (Int, Int) Int 993 nbc = strictArray (array (pHbounds dimension_all_zones) 994 ([(i, j)=: 1 | i <- [kmin],j <- [lmin + 1..lmax]]++ 995 [(i, j)=: 2 | i <- [kmax + 1],j <- [lmin + 1..lmax]]++ 996 [(i, j)=: 1 | j <- [lmin],i <- [kmin + 1..kmax]]++ 997 [(i, j)=: 1 | j <- [lmax + 1],i <- [lmin + 1..lmax]]++ 998 [(i, j)=: 4 | i <- [kmin],j <- [lmin]]++ 999 [(i, j)=: 4 | i <- [kmin],j <- [lmax + 1]]++ 1000 [(i, j)=: 4 | i <- [kmax + 1],j <- [lmin]]++ 1001 [(i, j)=: 4 | i <- [kmax + 1],j <- [lmax + 1]])) 1002 1003 pbb :: DblVector 1004 pbb = strictArray (array (1, 4) 1005 ([2 =: 6.0]++ 1006 [i=: 0.0 | i <- [1..4], not (i == 2)])) 1007 1008 pb:: DblVector 1009 pb = strictArray (array (1, 4) 1010 ([i=: 1.0 | i <- [1..4], not (i == 2 || i == 3)]++ 1011 [i=: 0.0 | i <- [2..3]])) 1012 1013 qb:: DblVector 1014 qb = pb 1015 1016 all_zero_nodes :: a -> DblArray 1017 all_zero_nodes _ = strictArray (array (pHbounds dimension_all_nodes) 1018 ([node=: 0.0 | node <- all_nodes])) 1019 1020 all_zero_zones :: a -> DblArray 1021 all_zero_zones _ = strictArray (array (pHbounds dimension_all_zones) 1022 ([zone=: 0.0 | zone <- all_zones])) 1023 1024 make_position_matrix :: ((Int, Int) -> (Double, Double)) -> (DblArray, DblArray) 1025 make_position_matrix interior_function = 1026 let { 1027 boundary_position rx zx ry zy ra za = 1028 let { 1029 (rax, zax) = (ra - rx, za - zx) ; 1030 (ryx, zyx) = (ry - rx, zy - zx) ; 1031 omega = ((2.0::Double) * (rax * ryx + zax * zyx)) 1032 / (ryx ^ (2::Int) + zyx ^ (2::Int)) ; 1033 rb = (rx - rax) + omega * ryx ; 1034 zb = (zx - zax) + omega * zyx 1035 } in (rb, zb) ; 1036 reflect_node x_dir y_dir a_dir node = 1037 let { 1038 rx = reflect x_dir node r' ; 1039 zx = reflect x_dir node z' ; 1040 ry = reflect y_dir node r' ; 1041 zy = reflect y_dir node z' ; 1042 ra = reflect a_dir node r' ; 1043 za = reflect a_dir node z' 1044 } in boundary_position rx zx ry zy ra 1045 za ; 1046 (r', z') = (arrays_2 (pHbounds dimension_all_nodes) 1047 ([node=: interior_function node 1048 | node <- interior_nodes]++ 1049 [node=: reflect_node south southeast 1050 farsouth node 1051 | node <- north_nodes]++ 1052 [node=: reflect_node north northeast 1053 farnorth node 1054 | node <- south_nodes]++ 1055 [node=: reflect_node west southwest 1056 farwest node | node <- east_nodes]++ 1057 [node=: reflect_node east southeast 1058 fareast node | node <- west_nodes]++ 1059 [node=: reflect_node southwest west 1060 farwest node | node <- north_east_corner]++ 1061 [node=: reflect_node northwest west 1062 farwest node | node <- south_east_corner]++ 1063 [node=: reflect_node southeast south 1064 farsouth node 1065 | node <- north_west_corner]++ 1066 [node=: reflect_node northeast east 1067 fareast node | node <- south_west_corner]++ 1068 [node=: reflect_node south southwest 1069 farsouth node 1070 | node <- west_of_north_east]++ 1071 [node=: reflect_node north northwest 1072 farnorth node 1073 | node <- west_of_south_east]++ 1074 [node=: reflect_node west northwest 1075 farwest node | node <- north_of_south_east]++ 1076 [node=: reflect_node east northeast 1077 fareast node | node <- north_of_south_west])) 1078 } in strictArrays_2 (r', z') 1079 1080 zone_area_vol :: (DblArray, DblArray) -> Zone -> (Double, Double) 1081 zone_area_vol (r, z) zone = 1082 let { 1083 (rsw, zsw) = (r!zone_corner_southwest zone, 1084 z!zone_corner_southwest zone) ; 1085 (rse, zse) = (r!zone_corner_southeast zone, 1086 z!zone_corner_southeast zone) ; 1087 (rne, zne) = (r!zone_corner_northeast zone, 1088 z!zone_corner_northeast zone) ; 1089 (rnw, znw) = (r!zone_corner_northwest zone, 1090 z!zone_corner_northwest zone) ; 1091 area1 = (rse * (zne - zsw) + rsw * (zse - zne)) 1092 + rne * (zsw - zse) ; 1093 area2 = (rnw * (zsw - zne) + rne * (znw - zsw)) 1094 + rsw * (zne - znw) ; 1095 average1 = (rsw + rse) + rne ; 1096 volume1 = (area1 * average1) / (3.0::Double) ; 1097 average2 = (rsw + rne) + rnw ; 1098 volume2 = (area2 * average2) / (3.0::Double) ; 1099 zone_area = (area1 + area2) / (2.0::Double) ; 1100 zone_volume = (volume1 + volume2) * (pi::Double) 1101 } in (zone_area, zone_volume) 1102 1103 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1104 -- %% Utilities 1105 -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1106 1107 reduce_list :: (a -> b -> a) -> a -> [b] -> a 1108 reduce_list fcn z list = foldl fcn z list 1109 1110 max_list:: [Double] -> Double 1111 max_list list = reduce_list max (head (list :: [ Double])) list 1112 1113 min_list:: [Double] -> Double 1114 min_list list = reduce_list min (head (list :: [ Double])) list 1115 1116 sum_list:: [Double] -> Double 1117 sum_list list = reduce_list (+) (0.0 :: Double) (list :: [ Double]) 1118 1119 -- %%%% Edit History: /home/prj2/rpaul/ph/tests/ek-simple/003/ek-simple.id 1120 1121 -- %%%+/tmp_mnt/home/prj2/rpaul/ph/tests/ek-simple-93.id 1122 -- %%% [rpaul] 4/13/94 17:28 verified computations: velocity, position, area/volume/density , artificial-viscosity and heat-conduction. 1123 -- %%% - fixed the area computation in zone_area. (equation 24) 1124 -- %%% - both equation 24 and the code are still missing the constant 1125 -- %%% factor of 2pi in the volume computation 1126 -- %%% - the code for artifical viscocity (equation 27) were missing upper_disc^2 1127 -- %%% and lower_disc^2 1128 -- %%% - the code for interior_cc (equation 37) was missing theta_hat^5/2 1129 -- %%% - the initial condition for heat is now 10.0 (was .0001)