Adapt `nofib` code to Foldable-generalised Prelude
[nofib.git] / real / fluid / Jcb_method.hs
1 {-
2 Implementation of the Jacobi iteration
3
4 XZ, 24/10/91
5 -}
6
7 {-
8 Modified to adopt S_array
9
10 XZ, 19/2/92
11 -}
12
13 module Jcb_method ( jcb_method ) where
14
15 import Defs
16 import S_Array -- not needed w/ proper module handling
17 import Norm -- ditto
18 import Asb_routs
19 import Tol_cal
20 import Data.Ix
21
22 -----------------------------------------------------------
23 -- Jacobi iteration method: --
24 -- solves: scalor*MX = B --
25 -- with iteration counter (r) as --
26 -- X(r+1) = X(r) + relax*(B/scalor - MX(r))/D --
27 -- The system matrix M is never really assembled --
28 -----------------------------------------------------------
29
30 jcb_method
31 :: Int
32 -> My_Array Int (Frac_type,((Frac_type,Frac_type,Frac_type),
33 (Frac_type,Frac_type,Frac_type)))
34 -> My_Array Int [(Int,Int)]
35 -> My_Array Int [Int]
36 -> (My_Array Int Bool,(My_Array Int Bool, My_Array Int Bool))
37 -> (My_Array Int Frac_type, My_Array Int Frac_type)
38 -> Frac_type
39 -> Int
40 -> Frac_type
41 -> Frac_type
42 -> (My_Array Int Frac_type, My_Array Int Frac_type)
43
44 sparse_elems = \y -> map (\(_,x)->x) (sparse_assocs y)
45
46 jcb_method f_step el_det_fac asb_table v_steer
47 (all_bry,(x_fixed,y_fixed)) (b1,b2) scalor
48 max_iter m_iter_toler relax =
49 sub_jacobi (init_u,init_u) max_iter
50 where
51 n_bnds = s_bounds asb_table
52 init_u = s_array n_bnds []
53 -- the recursive function of the Jacobi iteration
54 sub_jacobi old_x n = -- X(r) iteration_counter
55 if
56 -- checking the iteration limit
57 ( n <= 1 ) ||
58 -- checking the tolerance
59 (
60 (n /= max_iter) &&
61 ( tol_cal
62 ((s_elems (fst new_x)) ++ (s_elems (snd new_x)))
63 ((sparse_elems (fst diff)) ++ (sparse_elems (snd diff)))
64 True
65 ) < m_iter_toler
66 )
67 {-
68 (
69 (sum ( map abs ( s_elems (fst diff) ) )) +
70 (sum ( map abs ( s_elems (snd diff) ) ))
71 )
72 < m_iter_toler
73 -}
74 then new_x -- X(r+1)
75 else sub_jacobi new_x (n-1)
76 where
77 -- new adaptive x
78 new_x =
79 if ( n == max_iter )
80 -- first iteration
81 then diff
82 -- otherwise
83 else add_u old_x diff
84 diff =
85 if ( f_step == 1 )
86 then
87 -- For step 1. Only fixes some boundry nodes.
88 (
89 find_diff x_fixed (fst old_x) b1,
90 find_diff y_fixed (snd old_x) b2
91 )
92 else
93 -- For step 3. Fixes all boundry nodes.
94 (
95 find_diff all_bry (fst old_x) b1,
96 find_diff all_bry (snd old_x) b2
97 )
98 bindTo x k = k x -- old haskell 1.0 "let"
99
100 -- function which does one adaptation,
101 -- ie, calculates:
102 -- relax*(B - MX(r))/D
103 find_diff fixed x b = -- list_of_unfixed_nodes X(r) B
104 s_def_listArray n_bnds (0::Frac_type)
105 [
106 if fixed!^i
107 then 0
108 else
109 ((b!^i) / scalor) `bindTo` ( \ b' ->
110 ((
111 if ( n == max_iter )
112 then b'
113 else
114 b' -
115 sum [
116 (list_inner_prod
117 (get_val x (v_steer!^e))
118 (map (mult (fst (el_det_fac!^e))) (m_mat!^id)))
119 | (e,id)<-asb_table!^i
120 ]
121 ) * relax) / (mat_diag!^i) )
122 | i <- range n_bnds
123 ]
124 -- row sums of system M
125 mat_diag =
126 s_listArray n_bnds
127 [
128 sum
129 [ (fst ( el_det_fac!^e)) * (m_r_sum!^id)
130 | (e,id)<-asb_table!^i
131 ]
132 | i <- range n_bnds
133 ]
134
135 -----------------------------------------------------------
136 -- row sums of element matrix m_mat --
137 -----------------------------------------------------------
138
139 m_r_sum =
140 s_def_listArray (1,v_nodel) def_v
141 [ def_v,def_v,def_v,v',v',v' ]
142 where
143 def_v = (12 / 180)::Frac_type
144 v' = (68 / 180)::Frac_type
145
146 -----------------------------------------------------------
147 -- element matrix --
148 -----------------------------------------------------------
149
150 m_mat =
151 s_listArray (1,v_nodel)
152 (map (map (mult inv_180))
153 [
154 [6,(-1),(-1),(-4),0,0],
155 [(-1),6,(-1),0,(-4),0],
156 [(-1),(-1),6,0,0,(-4)],
157 [(-4),0,0,32,16,16],
158 [0,(-4),0,16,32,16],
159 [0,0,(-4),16,16,32]
160 ])
161 where
162 inv_180 = 1 / 180