1 {- 2 Taylor-Galerkin/Pressure-correction algorithm. 3 Solving four increamental matrix equations iteratively: 4 const1*M(U'-U) = rh1(U,P) 5 const2*M(U''-U) = rh1(U',P) 6 const3*K(P'-P) = rh2(U') 7 const4*M(U'''-U'') = rh3(P'-P) 8 The 3rd equation is solved by using the Choleski 9 decomposition menthod and the rest are solved by using 10 the Jacobi iteration method. 11 12 XZ, 24/10/91 13 -} 14 15 {- 16 Modified to adopt S_arrays. 17 18 Evaluations are forced by using normalize_obj. 19 (This is currently necessary for running the whole 20 program on a 200 element problem). 21 22 XZ, 19/2/92 23 -} 24 25 {- 26 Iteration along time-step implemented 27 28 XZ, 25/2/92 29 -} 30 31 module TG_iter ( tg_iter ) where 32 33 import Defs 34 import S_Array -- not needed w/ proper module handling 35 import Norm -- ditto 36 import Asb_routs 37 import Rhs_Asb_routs 38 import Jcb_method 39 import Chl_method 40 import Tol_cal 41 42 ----------------------------------------------------------- 43 -- Iterative TG algorithm. -- 44 -- Functions called : -- 45 -- for equation solving: "chl_method" and "jcb_method" -- 46 -- for RHS assembling : "get_rh1", "get_rh2" and -- 47 -- "get_rh3" -- 48 ----------------------------------------------------------- 49 50 tg_iter 51 :: Bool -> Int -> Float -> Int -> Float -> Float -> Float 52 -> (S_array (Float, ((Float, Float, Float), (Float, Float, Float)))) 53 -> (S_array [(Int, Int)]) 54 -> (S_array [(Int, Int)]) 55 -> (S_array [Int]) 56 -> (S_array [Int]) 57 -> (S_array Bool, (S_array Bool, S_array Bool)) 58 -> [Int] 59 -> (S_array (S_array Float, S_array (Int, [Float])), S_array Int) 60 -> (S_array Float, (S_array Float, S_array Float)) 61 -> String 62 63 tg_iter 64 mon m_iter m_toler max_jcb_iter jcb_toler 65 relax dlt_t el_det_fac v_asb_table p_asb_table 66 v_steer p_steer bry_nodes p_fixed chl_fac (p,u) = 67 do_tg_iter m_iter (pack_obj p,pack_obj u) [] 68 where 69 do_tg_iter n (old_p,old_u) res = 70 if (n<=1) || (max_tol<m_toler) 71 then new_res ++ show_final 72 else 73 if max_tol>large 74 then error "main iteration: overflow!" 75 else do_tg_iter (n-1) (new_p,new_u) new_res 76 where 77 new_res = 78 res ++ 79 "at time step " ++ (shows (m_iter-n+1) ":\n\n") ++ 80 (if mon 81 then 82 "initial presure:\n" ++ (shows (retrieve_obj old_p) "\n") ++ 83 "inititial velocities:\n" ++ (shows (retrieve_obj old_u) "\n") ++ 84 "rh1a:\n" ++ (shows rh1a "\n") ++ 85 "velocities after the 1st half of step 1:\n" ++ (shows (retrieve_obj u1a) "\n") ++ 86 "rh1b:\n" ++ (shows rh1b "\n") ++ 87 "velocities after the 2nd half of step 1:\n" ++ (shows (retrieve_obj u1b) "\n") ++ 88 "rh2:\n" ++ (shows (retrieve_obj rh2) "\n") ++ 89 "rh3:\n" ++ (shows rh3 "\n") 90 else "") 91 show_final = 92 "presure:\n" ++ (shows (retrieve_obj new_p) "\n") ++ 93 "velocities:\n" ++ (shows (retrieve_obj new_u) "\n") 94 tmp_p | normalize_obj new_p = retrieve_obj new_p 95 (tmp_u_x,tmp_u_y) | normalize_obj new_u = retrieve_obj new_u 96 max_tol = max tol_u tol_p 97 tol_u = 98 tol_cal ((s_elems tmp_u_x)++(s_elems tmp_u_y)) 99 ((subs tmp_u_x old_u_x) ++ 100 (subs tmp_u_y old_u_y)) 101 False 102 where 103 (old_u_x,old_u_y) | normalize_obj old_u = retrieve_obj old_u 104 tol_p | normalize_obj old_p = 105 tol_cal (s_elems tmp_p) 106 (subs tmp_p (retrieve_obj old_p)) False 107 -- step 1a 108 rh1a | normalize_obj old_p && normalize_obj old_u = 109 get_rh1 el_det_fac v_asb_table v_steer p_steer 110 (retrieve_obj old_p,retrieve_obj old_u) 111 del_u1a = 112 jcb_method 1 el_det_fac v_asb_table v_steer 113 bry_nodes rh1a (2/dlt_t) max_jcb_iter jcb_toler relax 114 u1a = pack_obj (add_u (retrieve_obj old_u) del_u1a) 115 -- step 1b 116 rh1b | normalize_obj u1a = 117 get_rh1 el_det_fac v_asb_table v_steer p_steer 118 (retrieve_obj old_p,retrieve_obj u1a) 119 del_u1b = 120 jcb_method 1 el_det_fac v_asb_table v_steer 121 bry_nodes rh1b (1/dlt_t) max_jcb_iter jcb_toler relax 122 u1b = pack_obj (add_u (retrieve_obj old_u) del_u1b) 123 -- step 2 124 rh2 | normalize_obj u1b = 125 pack_obj 126 (get_rh2 el_det_fac p_asb_table v_steer p_fixed 127 (retrieve_obj u1b)) 128 del_p | normalize_obj rh2 = 129 pack_obj (chl_method chl_fac (retrieve_obj rh2) (2/dlt_t)) 130 new_p = 131 pack_obj (add_mat (retrieve_obj old_p) (retrieve_obj del_p)) 132 -- step 3 133 rh3 | normalize_obj del_p = 134 get_rh3 el_det_fac v_asb_table p_steer (retrieve_obj del_p) 135 del_u3 = 136 jcb_method 3 el_det_fac v_asb_table v_steer 137 bry_nodes rh3 (2/dlt_t) max_jcb_iter jcb_toler relax 138 new_u = pack_obj (add_u (retrieve_obj u1b) del_u3) 139 140 subs = \x' x -> zipWith (-) (s_elems x') (s_elems x)