1 | from anuga.abstract_2d_finite_volumes.generic_domain import Generic_Domain |
---|
2 | import anuga.abstract_2d_finite_volumes.mesh_factory as mesh_factory |
---|
3 | from anuga.utilities.sparse import Sparse, Sparse_CSR |
---|
4 | from anuga.utilities.cg_solve import conjugate_gradient |
---|
5 | import anuga.abstract_2d_finite_volumes.neighbour_mesh as neighbour_mesh |
---|
6 | from generic_boundary_conditions import Dirichlet_boundary |
---|
7 | import numpy as num |
---|
8 | import anuga.utilities.log as log |
---|
9 | """ |
---|
10 | Kinematic Viscosity Class for ANUGA |
---|
11 | Lindon Roberts, 2010 |
---|
12 | |
---|
13 | Possible extensions: |
---|
14 | (1) Currently assume that the line between centroids |
---|
15 | is parallel to the outward normal. Change formula |
---|
16 | for this. |
---|
17 | (2) Higher order approximation (eg. Least squares |
---|
18 | fitting or quadratic etc.) |
---|
19 | (3) Vectorise the for loops - multiply by diagonal matrices? |
---|
20 | (4) Tidy up direct application calling structure: __mul__ for matrices, |
---|
21 | then ensure correct setting of include_boundaries |
---|
22 | (5) More comprehensive testing (larger problems) |
---|
23 | """ |
---|
24 | |
---|
25 | |
---|
26 | """ |
---|
27 | How to use: |
---|
28 | |
---|
29 | (1) Direct application |
---|
30 | domain = anuga.abstract_2d_finite_volumes.generic_domain.Generic_Domain(...) |
---|
31 | KV_Operator = Kinematic_Viscosity_Operator(domain) |
---|
32 | h = [n*1 numpy column vector] |
---|
33 | KV_Operator.apply_stage_heights(h) |
---|
34 | uh = [n*1 numpy column vector] |
---|
35 | vh = [n*1 numpy column vector] |
---|
36 | KV_Operator.set_qty_considered('u') #or KV_Operator.set_qty_considered(1) |
---|
37 | KV_Operator.build_boundary_vector() |
---|
38 | div_u = KV_Operator.apply_kv_operator(uh) |
---|
39 | KV_Operator.set_qty_considered('v') #or KV_Operator.set_qty_considered(2) |
---|
40 | KV_Operator.build_boundary_vector() |
---|
41 | div_v = KV_Operator.apply_kv_operator(vh) |
---|
42 | |
---|
43 | [Now div_u and div_v are (n,) numpy vectors] |
---|
44 | |
---|
45 | ------------------------------------------------------------------------ |
---|
46 | (2) Implicit use |
---|
47 | domain = anuga.abstract_2d_finite_volumes.generic_domain.Generic_Domain(...) |
---|
48 | KV_Operator = Kinematic_Viscosity_Operator(domain) |
---|
49 | h = [n*1 numpy column vector] |
---|
50 | KV_Operator.apply_stage_heights(h) |
---|
51 | div = [n*2 numpy matrix] #column 0 = div_u, column 1 = div_v |
---|
52 | velocities = KV_Operator.cg_solve(div) |
---|
53 | |
---|
54 | [Now velocities is a (n,2) numpy matrix = [(uh) (vh)]] |
---|
55 | |
---|
56 | If you just want to solve for (uh), say: |
---|
57 | KV_Operator.apply_stage_heights(h) |
---|
58 | div_u = [(n,) numpy vector] |
---|
59 | uh = KV_Operator.cg_solve(div_u, 'u') #or KV_Operator.cg_solve(div_u,1) |
---|
60 | """ |
---|
61 | |
---|
62 | class Kinematic_Viscosity_Operator: |
---|
63 | ## @brief Initialise a KV operator |
---|
64 | # @param domain The anuga.abstract_2d_finite_volumes.generic_domain.Generic_Domain for the problem |
---|
65 | # @param triangle_areas Whether we should scale the divergence by 1/|Triangle_i| |
---|
66 | # @param verbose |
---|
67 | def __init__(self, domain, triangle_areas=True, verbose=False): |
---|
68 | #Expose the domain attributes |
---|
69 | self.domain = domain |
---|
70 | self.mesh = domain.mesh |
---|
71 | self.boundary = domain.boundary |
---|
72 | self.boundary_enum = self.enumerate_boundary() |
---|
73 | self.n = len(self.domain) |
---|
74 | self.boundary_len = len(domain.boundary) |
---|
75 | self.qty_considered = 1 #1 or 2 (uh or vh respectively) |
---|
76 | self.verbose = verbose |
---|
77 | if self.verbose: log.critical('Kinematic Viscosity: Building structure') |
---|
78 | self.geo_structure = self.build_geo_structure() #list of Sparse matrices |
---|
79 | self.operator_matrix = Sparse(self.n,self.n+self.boundary_len) |
---|
80 | self.stage_heights = num.zeros((self.n,1),num.float) |
---|
81 | self.stage_heights_scaling = Sparse(self.n,self.n) |
---|
82 | self.apply_triangle_areas = triangle_areas |
---|
83 | self.boundary_vector = num.zeros((self.n,),num.float) |
---|
84 | |
---|
85 | self.triangle_areas = Sparse(self.n,self.n) |
---|
86 | for i in range(self.n): |
---|
87 | self.triangle_areas[i,i] = 1.0 / self.mesh.areas[i] |
---|
88 | if self.verbose: log.critical('Kinematic Viscosity: Initialisation Done') |
---|
89 | |
---|
90 | ## |
---|
91 | # Enumerate the boundary conditions in some way |
---|
92 | def enumerate_boundary(self): |
---|
93 | #Just enumerate by the dictionary.keys() list order |
---|
94 | #NOTE If this is changed, need to fix build_boundary_vector() method |
---|
95 | enumeration = {} |
---|
96 | keys = self.boundary.keys() |
---|
97 | for k,key in enumerate(keys): |
---|
98 | enumeration[key] = k |
---|
99 | return enumeration |
---|
100 | |
---|
101 | #Create the matrix structure for the geometric component of the operator |
---|
102 | # (constant across all applications of the operator) |
---|
103 | def build_geo_structure(self): |
---|
104 | A = [Sparse(self.n,self.n+self.boundary_len), Sparse(self.n,self.n+self.boundary_len), Sparse(self.n,self.n+self.boundary_len)] |
---|
105 | for i in range(self.n): |
---|
106 | this_centroid_coords = self.mesh.centroid_coordinates[i,:] |
---|
107 | for edge in range(3): |
---|
108 | j = self.mesh.neighbours[i,edge] |
---|
109 | edge_length = self.mesh.edgelengths[i,edge] |
---|
110 | dist = 0.0 |
---|
111 | index = j |
---|
112 | if j==-1: #Boundary edge |
---|
113 | boundary_coords = self.mesh.get_edge_midpoint_coordinate(i,edge) |
---|
114 | dist = num.linalg.norm(this_centroid_coords - boundary_coords) |
---|
115 | index = self.n + self.boundary_enum[(i,edge)] |
---|
116 | else: #Interior edge |
---|
117 | other_centroid_coords = self.mesh.centroid_coordinates[j,:] |
---|
118 | dist = num.linalg.norm(this_centroid_coords - other_centroid_coords) |
---|
119 | A[edge][i,i] = -edge_length / dist |
---|
120 | A[edge][i,index] = edge_length / dist |
---|
121 | return A |
---|
122 | |
---|
123 | def set_qty_considered(self,qty): |
---|
124 | if qty == 1 or qty == 'u': |
---|
125 | self.qty_considered = 1 |
---|
126 | elif qty == 2 or qty == 'v': |
---|
127 | self.qty_considered = 2 |
---|
128 | else: #Raise an exception |
---|
129 | msg = "Incorrect input qty" |
---|
130 | assert 0==1, msg |
---|
131 | |
---|
132 | def apply_stage_heights(self,h): |
---|
133 | msg = "(KV_Operator.apply_stage_heights) h vector has incorrect length" |
---|
134 | assert h.size == self.n, msg |
---|
135 | |
---|
136 | #h is a vector of self.n stage heights (1 per triangle) |
---|
137 | #scaled_geo_structure = num.zeros((3,self.n,self.n+self.boundary_len),num.float) |
---|
138 | scaled_geo_structure = [Sparse(self.n,self.n+self.boundary_len), Sparse(self.n,self.n+self.boundary_len), Sparse(self.n,self.n+self.boundary_len)] |
---|
139 | for i in range(self.n): |
---|
140 | h_i = h[i] |
---|
141 | for edge in range(3): |
---|
142 | j = self.mesh.neighbours[i,edge] |
---|
143 | h_j = 0.0 |
---|
144 | index = j |
---|
145 | if j==-1: #Boundary edge |
---|
146 | h_j = self.boundary[(i,edge)].evaluate()[0] |
---|
147 | index = self.n + self.boundary_enum[(i,edge)] |
---|
148 | else: #Interior edge |
---|
149 | h_j = h[j] |
---|
150 | scaled_geo_structure[edge][i,i] = 0.5*(h_i+h_j)*self.geo_structure[edge][i,i] |
---|
151 | scaled_geo_structure[edge][i,index] = 0.5*(h_i+h_j)*self.geo_structure[edge][i,index] |
---|
152 | self.operator_matrix = Sparse_CSR(scaled_geo_structure[0]+scaled_geo_structure[1]+scaled_geo_structure[2]) |
---|
153 | #self.operator_matrix = scaled_geo_structure[0]+scaled_geo_structure[1]+scaled_geo_structure[2] |
---|
154 | self.stage_heights = h |
---|
155 | #Set up the scaling matrix |
---|
156 | for i in range(self.n): |
---|
157 | if h[i] == 0.0: |
---|
158 | self.stage_heights_scaling[i,i] = 0.0 |
---|
159 | else: |
---|
160 | self.stage_heights_scaling[i,i] = 1.0/h[i] |
---|
161 | return self.operator_matrix |
---|
162 | |
---|
163 | def build_boundary_vector(self): |
---|
164 | #Require self.qty_considered already set and stage heights applied |
---|
165 | X = num.zeros((self.n+self.boundary_len,),num.float) |
---|
166 | #TODO Can vectorise? |
---|
167 | for i,key in enumerate(self.boundary.keys()): |
---|
168 | quantities = self.boundary[key].evaluate() |
---|
169 | h_i = quantities[0] |
---|
170 | if h_i == 0.0: |
---|
171 | X[self.n+i] = 0.0 |
---|
172 | else: |
---|
173 | X[self.n+i] = quantities[self.qty_considered] / h_i |
---|
174 | self.boundary_vector = self.operator_matrix * X |
---|
175 | #Tidy up |
---|
176 | if self.apply_triangle_areas: |
---|
177 | self.boundary_vector = self.triangle_areas * self.boundary_vector |
---|
178 | |
---|
179 | return self.boundary_vector |
---|
180 | |
---|
181 | def apply_kv_operator(self,V,qty_considered=None,include_boundary=True): |
---|
182 | msg = "(KV_Operator.apply_vector) V vector has incorrect dimensions" |
---|
183 | assert V.shape == (self.n,) or V.shape == (self.n,1), msg |
---|
184 | |
---|
185 | if not qty_considered==None: |
---|
186 | print "Quantity Considered changed!" |
---|
187 | self.set_qty_considered(qty_considered) |
---|
188 | self.build_boundary_vector() |
---|
189 | |
---|
190 | D = num.zeros((self.n,),num.float) |
---|
191 | #Change (uh) to (u), (vh) to (v) |
---|
192 | X = self.stage_heights_scaling * V |
---|
193 | |
---|
194 | if self.apply_triangle_areas: |
---|
195 | X = self.triangle_areas * X |
---|
196 | |
---|
197 | if include_boundary: |
---|
198 | #TODO Tidy this up |
---|
199 | #A = self.operator_matrix.todense() |
---|
200 | #D = num.mat(A[:,:self.n]) * num.mat(X.reshape(self.n,1)) + self.boundary_vector |
---|
201 | A = self.operator_matrix |
---|
202 | D = A * num.concatenate((X, num.zeros(self.boundary_len,num.float))) + self.boundary_vector |
---|
203 | else: |
---|
204 | A = self.operator_matrix |
---|
205 | D = A * num.concatenate((X, num.zeros(self.boundary_len,num.float))) |
---|
206 | #A = self.operator_matrix.todense() |
---|
207 | #D = num.mat(A[:,:self.n]) * num.mat(X.reshape(self.n,1)) |
---|
208 | |
---|
209 | return num.array(D).reshape(self.n,) |
---|
210 | |
---|
211 | def __mul__(self,other): |
---|
212 | try: |
---|
213 | B = num.array(other) |
---|
214 | except: |
---|
215 | msg = "Trying to multiply the Kinematic Viscosity Operator against a non-numeric type" |
---|
216 | raise msg |
---|
217 | |
---|
218 | if len(B.shape) == 0: |
---|
219 | #Scalar |
---|
220 | R = B*self |
---|
221 | elif len(B.shape) == 1 or (len(B.shape) == 2 and B.shape[1]==1): |
---|
222 | #Vector |
---|
223 | #include_boundary=False is this is *only* used for cg_solve() |
---|
224 | R = self.apply_kv_operator(other,include_boundary=False) |
---|
225 | else: |
---|
226 | raise ValueError, 'Dimension too high: d=%d' %len(B.shape) |
---|
227 | return R |
---|
228 | |
---|
229 | def __rmul__(self, other): |
---|
230 | #Right multiply with scalar |
---|
231 | try: |
---|
232 | other = float(other) |
---|
233 | except: |
---|
234 | msg = 'Sparse matrix can only "right-multiply" onto a scalar' |
---|
235 | raise TypeError, msg |
---|
236 | else: |
---|
237 | new = self.operator_matrix * new |
---|
238 | return new |
---|
239 | |
---|
240 | def cg_solve(self,B,qty_to_consider=None): |
---|
241 | if len(B.shape) == 1: |
---|
242 | return self.cg_solve_vector(B,qty_to_consider) |
---|
243 | elif len(B.shape) == 2: |
---|
244 | return self.cg_solve_matrix(B) |
---|
245 | else: |
---|
246 | raise ValueError, 'Dimension too high: d=%d' %len(B.shape) |
---|
247 | |
---|
248 | def cg_solve_matrix(self,B): |
---|
249 | assert B.shape[1] < 3, "Input matrix has too many columns (max 2)" |
---|
250 | |
---|
251 | X = num.zeros(B.shape,num.float) |
---|
252 | for i in range(B.shape[1]): |
---|
253 | X[:,i] = self.cg_solve_vector(B[:,i],i+1) #assuming B columns are (uh) and (vh) |
---|
254 | return X |
---|
255 | |
---|
256 | def cg_solve_vector(self,b,qty_to_consider=None): |
---|
257 | if not qty_to_consider==None: |
---|
258 | self.set_qty_considered(qty_to_consider) |
---|
259 | self.build_boundary_vector() |
---|
260 | if self.verbose : |
---|
261 | iprint = 1 |
---|
262 | else: |
---|
263 | iprint = 0 |
---|
264 | x = conjugate_gradient(self,b - self.boundary_vector, iprint=iprint) #Call the ANUGA conjugate gradient utility |
---|
265 | return x |
---|
266 | |
---|
267 | ################################################################################ |
---|
268 | |
---|
269 | if __name__ == "__main__": |
---|
270 | #Set up a rectangular grid and check minimum/maximum principle |
---|
271 | |
---|
272 | import time |
---|
273 | t0 = time.time() |
---|
274 | points,elements,boundary = mesh_factory.rectangular_cross(10, 15) |
---|
275 | mesh = neighbour_mesh.Mesh(points,elements) |
---|
276 | n = len(mesh) |
---|
277 | boundary_map = {} |
---|
278 | for (vol_id,edge) in boundary.keys(): |
---|
279 | x = mesh.get_edge_midpoint_coordinate(vol_id,edge)[0] |
---|
280 | y = mesh.get_edge_midpoint_coordinate(vol_id,edge)[1] |
---|
281 | #Define boundary conditions: u=x^2-y^2, v=x-y+2 |
---|
282 | boundary_map[(vol_id,edge)] = Dirichlet_boundary([1,x*x-y*y,x-y+2]) |
---|
283 | #Now define an operator |
---|
284 | domain = Generic_Domain(source=points,triangles=elements,boundary=boundary_map) |
---|
285 | KV_Operator = Kinematic_Viscosity_Operator(domain, verbose = True) |
---|
286 | |
---|
287 | #Set up the operator and RHS |
---|
288 | h = num.ones((n,),num.float) |
---|
289 | KV_Operator.apply_stage_heights(h) |
---|
290 | rhs = num.zeros((n,2),num.float) #we will solve for u and v |
---|
291 | |
---|
292 | #Solve |
---|
293 | velocities = KV_Operator.cg_solve(rhs) |
---|
294 | #Alternatively, velocities = [u v] where: |
---|
295 | u = KV_Operator.cg_solve(num.zeros((n,),num.float),'u') |
---|
296 | v = KV_Operator.cg_solve(num.zeros((n,),num.float),'v') |
---|
297 | vel_min = velocities.min(axis=0) |
---|
298 | vel_max = velocities.max(axis=0) |
---|
299 | |
---|
300 | #Check a minimum/maximum principle |
---|
301 | boundary_velocities = num.zeros((len(boundary_map),2),num.float) |
---|
302 | for i,key in enumerate(boundary_map.keys()): |
---|
303 | quantity = boundary_map[key].evaluate() |
---|
304 | boundary_velocities[i,:] = quantity[1:] |
---|
305 | boundary_vel_min = boundary_velocities.min(axis=0) |
---|
306 | boundary_vel_max = boundary_velocities.max(axis=0) |
---|
307 | |
---|
308 | |
---|
309 | print "Time ",time.time() - t0 |
---|
310 | print "Kinematic Viscosity Min/Max Principle Test" |
---|
311 | print "Minimum Principle?", (vel_min > boundary_vel_min).all() |
---|
312 | print "Maximum Principle?", (boundary_vel_max > vel_max).all() |
---|