source: trunk/anuga_work/development/2010-projects/kv/kv_solver_new.py @ 7884

Last change on this file since 7884 was 7884, checked in by steve, 14 years ago

Moving 2010 project

  • Property svn:executable set to *
File size: 12.7 KB
Line 
1from anuga.abstract_2d_finite_volumes.generic_domain import Generic_Domain
2import anuga.abstract_2d_finite_volumes.mesh_factory as mesh_factory
3from anuga.utilities.sparse import Sparse, Sparse_CSR
4from anuga.utilities.cg_solve import conjugate_gradient
5import anuga.abstract_2d_finite_volumes.neighbour_mesh as neighbour_mesh
6from generic_boundary_conditions import Dirichlet_boundary
7import numpy as num
8import anuga.utilities.log as log
9"""
10Kinematic Viscosity Class for ANUGA
11Lindon Roberts, 2010
12
13Possible extensions:
14(1) Currently assume that the line between centroids
15is parallel to the outward normal. Change formula
16for this.
17(2) Higher order approximation (eg. Least squares
18fitting or quadratic etc.)
19(3) Vectorise the for loops - multiply by diagonal matrices?
20(4) Tidy up direct application calling structure: __mul__ for matrices,
21then ensure correct setting of include_boundaries
22(5) More comprehensive testing (larger problems)
23"""
24
25
26"""
27How to use:
28
29(1) Direct application
30domain = anuga.abstract_2d_finite_volumes.generic_domain.Generic_Domain(...)
31KV_Operator = Kinematic_Viscosity_Operator(domain)
32h = [n*1 numpy column vector]
33KV_Operator.apply_stage_heights(h)
34uh = [n*1 numpy column vector]
35vh = [n*1 numpy column vector]
36KV_Operator.set_qty_considered('u') #or KV_Operator.set_qty_considered(1)
37KV_Operator.build_boundary_vector()
38div_u = KV_Operator.apply_kv_operator(uh)
39KV_Operator.set_qty_considered('v') #or KV_Operator.set_qty_considered(2)
40KV_Operator.build_boundary_vector()
41div_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
47domain = anuga.abstract_2d_finite_volumes.generic_domain.Generic_Domain(...)
48KV_Operator = Kinematic_Viscosity_Operator(domain)
49h = [n*1 numpy column vector]
50KV_Operator.apply_stage_heights(h)
51div = [n*2 numpy matrix] #column 0 = div_u, column 1 = div_v
52velocities = KV_Operator.cg_solve(div)
53
54[Now velocities is a (n,2) numpy matrix = [(uh) (vh)]]
55
56If you just want to solve for (uh), say:
57KV_Operator.apply_stage_heights(h)
58div_u = [(n,) numpy vector]
59uh = KV_Operator.cg_solve(div_u, 'u') #or KV_Operator.cg_solve(div_u,1)
60"""
61
62class 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
269if __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()
Note: See TracBrowser for help on using the repository browser.