# 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, 13 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
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
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.