Changeset 8152
- Timestamp:
- Mar 15, 2011, 6:02:19 PM (14 years ago)
- Location:
- trunk/anuga_core/source/anuga/operators
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/operators/kinematic_viscosity.py
r8135 r8152 1 2 1 from anuga import Domain 3 2 from anuga import Quantity … … 36 35 """ 37 36 38 def __init__(self, domain, diffusivity = None, triangle_areas=True, verbose=False):37 def __init__(self, domain, diffusivity = None, use_triangle_areas=True, verbose=False): 39 38 if verbose: log.critical('Kinematic Viscosity: Beginning Initialisation') 40 39 #Expose the domain attributes … … 43 42 self.mesh = domain.mesh 44 43 self.boundary = domain.boundary 45 44 self.boundary_enumeration = domain.boundary_enumeration 45 46 46 # Setup diffusivity quantity 47 47 if diffusivity is None: 48 48 diffusivity = Quantity(domain) 49 49 diffusivity.set_values(1.0) 50 diffusivity.set_boundary_values(1.0) 50 51 #check that diffusivity is a quantity associated with domain 51 52 assert diffusivity.domain == domain 52 53 53 54 self.diffusivity = diffusivity 54 self.diffusivity_bdry_data = self.diffusivity.boundary_values55 self.diffusivity_data = self.diffusivity.centroid_values55 #self.diffusivity_bdry_data = self.diffusivity.boundary_values 56 #self.diffusivity_cell_data = self.diffusivity.centroid_values 56 57 57 58 … … 61 62 self.tot_len = self.n + self.boundary_len 62 63 63 # FIXME SR: maybe this should already exist in mesh!64 self.boundary_enum = self.enumerate_boundary()65 66 64 self.verbose = verbose 67 65 … … 73 71 74 72 # Only needs to built once, doesn't change 75 kinematic_viscosity_ext.build_geo_structure(self , self.n, self.tot_len)73 kinematic_viscosity_ext.build_geo_structure(self) 76 74 77 75 # Setup type of scaling 78 self. apply_triangle_areas = triangle_areas76 self.set_triangle_areas(use_triangle_areas) 79 77 80 78 # FIXME SR: should this really be a matrix? … … 84 82 85 83 self.triangle_areas = Sparse_CSR(temp) 84 #self.triangle_areas 86 85 87 86 # FIXME SR: More to do with solving equation 88 87 self.qty_considered = 1 #1 or 2 (uh or vh respectively) 89 90 #FIXME SR: Do we need Sparse or should we just go with Sparse_CSR91 self.operator_matrix = Sparse(self.n, self.tot_len)92 88 93 89 #Sparse_CSR.data … … 98 94 self.operator_rowptr = 4 * num.arange(self.n + 1) 99 95 100 # Build matrix self.elliptic_ operator_matrix [A B]101 self.build_elliptic_ operator()96 # Build matrix self.elliptic_matrix [A B] 97 self.build_elliptic_matrix() 102 98 103 99 # Build self.boundary_term 104 self.build_operator_boundary_term()100 #self.build_elliptic_boundary_term() 105 101 106 102 self.parabolic_solve = False #Are we doing a parabolic solve at the moment? … … 108 104 if verbose: log.critical('Kinematic Viscosity: Initialisation Done') 109 105 106 def set_triangle_areas(self,flag=True): 107 108 self.apply_triangle_areas = flag 110 109 111 def enumerate_boundary(self):112 #Enumerate by boundary index113 # FIXME SR: Probably should be part of mesh114 enumeration = {}115 for i in range(self.n):116 for edge in range(3):117 j = self.mesh.neighbours[i, edge]118 if j < 0:119 enumeration[(i, edge)] = -j-1120 return enumeration121 110 122 111 def set_qty_considered(self, qty): … … 130 119 assert 0 == 1, msg 131 120 132 def build_elliptic_ operator(self):121 def build_elliptic_matrix(self): 133 122 """ 134 123 Builds matrix representing … … 141 130 #Arrays self.operator_data, self.operator_colind, self.operator_rowptr 142 131 # are setup via this call 143 kinematic_viscosity_ext.build_operator_matrix(self,self.diffusivity_data, self.diffusivity_bdry_data) 144 145 self.elliptic_operator_matrix = Sparse_CSR(None, \ 132 kinematic_viscosity_ext.build_elliptic_matrix(self, \ 133 self.diffusivity.centroid_values, \ 134 self.diffusivity.boundary_values) 135 136 self.elliptic_matrix = Sparse_CSR(None, \ 146 137 self.operator_data, self.operator_colind, self.operator_rowptr, \ 147 138 self.n, self.tot_len) 139 140 #print 'elliptic_matrix' 141 #print self.elliptic_matrix 148 142 149 143 # #Set up the scaling matrix … … 153 147 # Sparse_CSR(None, self.diffusivity_data, num.arange(self.n), num.arange(self.n +1), self.n, self.n) 154 148 155 def update_elliptic_ operator(self):149 def update_elliptic_matrix(self): 156 150 """ 157 151 Updates the data values of matrix representing … … 165 159 # through to the Sparse_CSR matrix. 166 160 167 kinematic_viscosity_ext.update_operator_matrix(self, self.n, self.tot_len, \ 168 self.diffusivity_data, self.diffusivity_bdry_data) 169 170 171 172 173 def build_operator_boundary_term(self): 161 kinematic_viscosity_ext.update_elliptic_matrix(self, \ 162 self.diffusivity.centroid_values, \ 163 self.diffusivity.boundary_values) 164 165 166 def build_elliptic_boundary_term(self,quantity): 174 167 """ 175 168 Operator has form [A B] and U = [ u ; b] … … 187 180 X = num.zeros((tot_len,), num.float) 188 181 189 X[n:] = self.diffusivity_bdry_data190 self.boundary_term[:] = self. operator_matrix * X182 X[n:] = quantity.boundary_values 183 self.boundary_term[:] = self.elliptic_matrix * X 191 184 192 185 #Tidy up … … 195 188 196 189 197 def elliptic_multiply(self, V, qty_considered=None, include_boundary=True): 198 msg = "(KV_Operator.apply_vector) V vector has incorrect dimensions" 199 assert V.shape == (self.n, ) or V.shape == (self.n, 1), msg 200 201 if qty_considered != None: 202 print "Quantity Considered changed!" 203 self.set_qty_considered(qty_considered) 204 205 D = num.zeros((self.n, ), num.float) 206 #Change (uh) to (u), (vh) to (v) 207 X = self.stage_heights_scaling * V 190 def elliptic_multiply(self, quantity_in, quantity_out=None, include_boundary=True): 191 192 n = self.n 193 tot_len = self.tot_len 194 195 V = num.zeros((tot_len,), num.float) 196 X = num.zeros((n,), num.float) 197 198 if quantity_out is None: 199 quantity_out = Quantity(self.domain) 200 201 V[0:n] = quantity_in.centroid_values 202 V[n:] = 0.0 203 204 205 # FIXME SR: These sparse matrix vector multiplications 206 # should be done in such a way to reuse memory, ie 207 # should have a procedure 208 # matrix.multiply(vector_in, vector_out) 209 208 210 209 211 if self.apply_triangle_areas: 210 X = self.triangle_areas * X211 212 X = num.append(X, num.zeros((self.boundary_len, )), axis=0) 213 D = self.operator_matrix * X212 V[0:n] = self.triangle_areas * V[0:n] 213 214 215 X[:] = self.elliptic_matrix * V 214 216 215 217 if include_boundary: 216 D += self.boundary_vector[:, self.qty_considered-1] 217 #make sure we return (uh) not (u), for the solver's benefit 218 D = self.stage_heights * D 219 220 return num.array(D).reshape(self.n, ) 218 self.build_elliptic_boundary_term(quantity_in) 219 X[:] += self.boundary_term 220 221 222 quantity_out.set_values(X, location = 'centroids') 223 return quantity_out 221 224 222 225 #parabolic_multiply(V) = identity - dt*elliptic_multiply(V) … … 237 240 #Multiply out 238 241 X = num.append(X, num.zeros((self.boundary_len, )), axis=0) 239 D = D - self.dt * (self. operator_matrix * X)242 D = D - self.dt * (self.elliptic_matrix * X) 240 243 241 244 return num.array(D).reshape(self.n, ) … … 270 273 raise TypeError, msg 271 274 else: 272 new = self. operator_matrix * new275 new = self.elliptic_matrix * new 273 276 return new 274 277 -
trunk/anuga_core/source/anuga/operators/kinematic_viscosity_ext.c
r8135 r8152 63 63 // Edge 64 64 if (j < 0 ) { 65 65 m = -j - 1; 66 66 geo_indices[3*i+edge] = n + m; 67 67 edge_counted++; … … 86 86 87 87 88 int build_ operator_matrix(int n,88 int build_elliptic_matrix(int n, 89 89 int tot_len, 90 90 long *geo_indices, … … 132 132 } 133 133 134 int update_ operator_matrix(int n,134 int update_elliptic_matrix(int n, 135 135 int tot_len, 136 136 long *geo_indices, … … 194 194 195 195 //Convert Python arguments to C 196 if (!PyArg_ParseTuple(args, "O ii", &kv_operator)) {196 if (!PyArg_ParseTuple(args, "O", &kv_operator)) { 197 197 PyErr_SetString(PyExc_RuntimeError, "build_geo_structure could not parse input"); 198 198 return NULL; … … 241 241 } 242 242 243 static PyObject *py_build_ operator_matrix(PyObject *self, PyObject *args) {243 static PyObject *py_build_elliptic_matrix(PyObject *self, PyObject *args) { 244 244 PyObject *kv_operator; 245 245 int n, tot_len, err; … … 254 254 //Convert Python arguments to C 255 255 if (!PyArg_ParseTuple(args, "OOO", &kv_operator, &cell_data, &bdry_data)) { 256 PyErr_SetString(PyExc_RuntimeError, "build_ operator_matrix could not parse input");256 PyErr_SetString(PyExc_RuntimeError, "build_elliptic_matrix could not parse input"); 257 257 return NULL; 258 258 } … … 268 268 colind = get_consecutive_array(kv_operator,"operator_colind"); 269 269 270 err = build_ operator_matrix(n,tot_len,270 err = build_elliptic_matrix(n,tot_len, 271 271 (long *)geo_indices -> data, 272 272 (double *)geo_values -> data, … … 289 289 290 290 291 static PyObject *py_update_ operator_matrix(PyObject *self, PyObject *args) {291 static PyObject *py_update_elliptic_matrix(PyObject *self, PyObject *args) { 292 292 PyObject *kv_operator; 293 293 int n, tot_len, err; … … 302 302 //Convert Python arguments to C 303 303 if (!PyArg_ParseTuple(args, "OOO", &kv_operator, &cell_data, &bdry_data)) { 304 PyErr_SetString(PyExc_RuntimeError, "update_ operator_matrix could not parse input");304 PyErr_SetString(PyExc_RuntimeError, "update_elliptic_matrix could not parse input"); 305 305 return NULL; 306 306 } … … 316 316 colind = get_consecutive_array(kv_operator,"operator_colind"); 317 317 318 err = update_ operator_matrix(n,tot_len,318 err = update_elliptic_matrix(n,tot_len, 319 319 (long *)geo_indices -> data, 320 320 (double *)geo_values -> data, … … 339 339 static struct PyMethodDef MethodTable[] = { 340 340 {"build_geo_structure",py_build_geo_structure,METH_VARARGS,"Print out"}, 341 {"build_ operator_matrix",py_build_operator_matrix,METH_VARARGS,"Print out"},342 {"update_ operator_matrix",py_update_operator_matrix,METH_VARARGS,"Print out"},341 {"build_elliptic_matrix",py_build_elliptic_matrix,METH_VARARGS,"Print out"}, 342 {"update_elliptic_matrix",py_update_elliptic_matrix,METH_VARARGS,"Print out"}, 343 343 {NULL,NULL,0,NULL} // sentinel 344 344 }; -
trunk/anuga_core/source/anuga/operators/test_kinematic_viscosity.py
r8135 r8152 1 1 from anuga import Domain 2 from anuga import Quantity 2 3 from anuga import Dirichlet_boundary 3 4 from kinematic_viscosity import Kinematic_Viscosity … … 33 34 34 35 36 35 37 #print domain.quantities['stage'].vertex_values 36 38 … … 63 65 domain.update_boundary() 64 66 67 68 65 69 return Kinematic_Viscosity(domain) 66 70 67 71 def test_enumerate_boundary(self): 68 72 operator1 = self.operator1() 69 boundary_enum = operator1.boundary_enum70 71 assert boundary_enum [(0,0)] == 072 assert boundary_enum [(0,1)] == 173 assert boundary_enum [(0,2)] == 273 boundary_enumeration = operator1.domain.boundary_enumeration 74 75 assert boundary_enumeration[(0,0)] == 0 76 assert boundary_enumeration[(0,1)] == 1 77 assert boundary_enumeration[(0,2)] == 2 74 78 75 79 operator2 = self.operator2() 76 boundary_enum = operator2.boundary_enum77 78 79 assert boundary_enum [(0,1)] == 080 assert boundary_enum [(0,2)] == 181 assert boundary_enum [(1,0)] == 282 assert boundary_enum [(1,2)] == 380 boundary_enumeration = operator2.domain.boundary_enumeration 81 82 83 assert boundary_enumeration[(0,1)] == 0 84 assert boundary_enumeration[(0,2)] == 1 85 assert boundary_enumeration[(1,0)] == 2 86 assert boundary_enumeration[(1,2)] == 3 83 87 84 88 def test_geo_structure(self): … … 96 100 assert num.allclose(values, num.array([[-3.0,-6.0/sqrt(5),-6.0/sqrt(5)],[-6.0/sqrt(5),-3.0,-6.0/sqrt(5)]])) 97 101 98 def test_elliptic_matrix(self): 99 operator1 = self.operator1() 100 operator2 = self.operator2() 101 102 domain1 = operator1.domain 103 domain2 = operator2.domain 104 105 106 A = operator1.elliptic_operator_matrix 107 108 print A 109 print num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)]) 110 102 def test_elliptic_matrix_one_triangle(self): 103 104 operator = self.operator1() 105 domain = operator.domain 106 107 operator.update_elliptic_matrix() 108 109 A = operator.elliptic_matrix 110 111 111 assert num.allclose(A.todense(), num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)])) 112 112 113 114 A = operator1.elliptic_operator_matrix 115 assert num.allclose(A.todense(), 1.5*num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)])) 116 117 h = num.array([1.0, 1.0]) 118 A = operator2.apply_stage_heights(h) 119 #From test_kv_system_geometry 113 diffusivity = operator.diffusivity 114 diffusivity.set_values(10.0) 115 diffusivity.set_boundary_values(10.0) 116 117 operator.update_elliptic_matrix() 118 119 assert num.allclose(A.todense(), 10*num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)])) 120 121 122 def test_elliptic_matrix_two_triangles(self): 123 124 125 operator = self.operator2() 126 127 domain = operator.domain 128 diffusivity = operator.diffusivity 129 130 A = operator.elliptic_matrix 131 132 diffusivity.set_values(1.0) 133 diffusivity.set_boundary_values(1.0) 134 operator.update_elliptic_matrix() 135 120 136 A0 = num.array([[-3.0,3.0,0.0,0.0,0.0,0.0], 121 137 [0.0,-6.0/sqrt(5.0),0.0,0.0,6.0/sqrt(5.0),0.0]]) … … 124 140 A2 = num.array([[-6.0/sqrt(5.0),0.0,0.0,6.0/sqrt(5.0),0.0,0.0],\ 125 141 [0.0, -6.0/sqrt(5.0), 0.0, 0.0, 0.0, 6.0/sqrt(5.0)]]) 142 143 126 144 assert num.allclose(A.todense(), A0+A1+A2) 127 145 128 h = num.array([2.0, 1.0]) 129 A = operator2.apply_stage_heights(h) 146 diffusivity.set_values([2.0, 1.0], location = 'centroids') 147 diffusivity.set_boundary_values(1.0) 148 operator.update_elliptic_matrix() 149 150 A = operator.elliptic_matrix 151 152 130 153 assert num.allclose(A.todense()[0,:], 1.5*A0[0,:]+1.5*A1[0,:]+1.5*A2[0,:]) 131 154 assert num.allclose(A.todense()[1,:], A0[1,:]+1.5*A1[1,:]+A2[1,:]) 132 155 133 h = num.array([-2.0, -2.0]) 134 A = operator2.apply_stage_heights(h) 156 diffusivity.set_values([-2.0, -2.0], location = 'centroids') 157 diffusivity.set_boundary_values(1.0) 158 operator.update_elliptic_matrix() 159 135 160 assert num.allclose(A.todense()[0,:], -2*A0[0,:]-0.5*A1[0,:]-0.5*A2[0,:]) 136 161 assert num.allclose(A.todense()[1,:], -0.5*A0[1,:]-2*A1[1,:]-0.5*A2[1,:]) 137 162 138 def test_elliptic_multiply(self): 139 operator1 = self.operator1() 140 operator1.apply_stage_heights(num.array([[1.0]])) #h=1 141 operator1.build_boundary_vector() 142 V1 = num.array([2.0]) #(uh)=2 <- Centriod value 143 V2 = num.array([2.0]) #(vh)=2 <- Centroid value 144 A = num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)]) 145 U1 = num.array([[2.0],[2.0],[1.0],[1.0]]) #(uh) for centroid, 3 edges 146 U2 = num.array([[2.0],[1.0],[2.0],[0.0]]) #(vh) for centroid, 3 edge 147 148 operator1.set_qty_considered('u') 149 D1 = operator1.elliptic_multiply(V1) 150 assert num.allclose(D1, 2*num.array(num.mat(A)*num.mat(U1)).reshape(1,)) #2* for triangle_areas 151 152 operator1.set_qty_considered(2) 153 D2 = operator1.elliptic_multiply(V2) 154 assert num.allclose(D2, 2*num.array(num.mat(A)*num.mat(U2)).reshape(1,)) #2* for triangle_areas 163 def test_elliptic_multiply_include_boundary_one_triangle(self): 164 operator = self.operator1() 165 operator.set_triangle_areas(False) 166 167 print operator.apply_triangle_areas 168 169 n = operator.n 170 171 q_in = Quantity(operator.domain) 172 q_in.set_values(1.0) 173 q_in.set_boundary_values(1.0) 174 operator.update_elliptic_matrix() 175 176 177 A = num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)]) 178 179 180 q_1 = operator.elliptic_multiply(q_in) 181 182 q_2 = operator.elliptic_multiply(q_in, quantity_out = q_in) 183 184 assert id(q_in) == id(q_2) 185 186 assert num.allclose(q_1.centroid_values,q_2.centroid_values) 187 188 assert num.allclose( num.zeros((n,), num.float), q_1.centroid_values ) 189 190 #Now have different boundary values 191 192 q_in.set_values(1.0) 193 q_in.set_boundary_values(0.0) 194 operator.update_elliptic_matrix() 195 196 197 A = num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)]) 198 199 200 q_1 = operator.elliptic_multiply(q_in) 201 202 assert num.allclose( [-6.0-12.0/sqrt(5)], q_1.centroid_values ) 203 204 def test_elliptic_multiply_exclude_boundary_one_triangle(self): 205 operator = self.operator1() 206 operator.set_triangle_areas(False) 207 208 print operator.apply_triangle_areas 209 #n = operator.n 210 211 q_in = Quantity(operator.domain) 212 q_in.set_values(1.0) 213 q_in.set_boundary_values(1.0) 214 operator.update_elliptic_matrix() 215 216 217 A = num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)]) 218 219 220 q_1 = operator.elliptic_multiply(q_in, include_boundary=False) 221 222 223 assert num.allclose( [-6.0-12.0/sqrt(5)], q_1.centroid_values ) 224 225 def test_elliptic_multiply_include_boundary_one_triangle(self): 226 operator = self.operator1() 227 operator.set_triangle_areas(True) 228 229 n = operator.n 230 231 q_in = Quantity(operator.domain) 232 q_in.set_values(1.0) 233 q_in.set_boundary_values(1.0) 234 operator.update_elliptic_matrix() 235 236 237 A = num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)]) 238 239 240 q_1 = operator.elliptic_multiply(q_in) 241 242 q_2 = operator.elliptic_multiply(q_in, quantity_out = q_in) 243 244 assert id(q_in) == id(q_2) 245 246 assert num.allclose(q_1.centroid_values,q_2.centroid_values) 247 248 assert num.allclose( num.zeros((n,), num.float), q_1.centroid_values ) 249 250 #Now have different boundary values 251 252 q_in.set_values(1.0) 253 q_in.set_boundary_values(0.0) 254 operator.update_elliptic_matrix() 255 256 257 A = num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)]) 258 259 260 q_1 = operator.elliptic_multiply(q_in) 261 262 assert num.allclose( [-12.0-24.0/sqrt(5)], q_1.centroid_values ) 263 264 def test_elliptic_multiply_exclude_boundary_one_triangle(self): 265 operator = self.operator1() 266 operator.set_triangle_areas(True) 267 268 q_in = Quantity(operator.domain) 269 q_in.set_values(1.0) 270 q_in.set_boundary_values(1.0) 271 operator.update_elliptic_matrix() 272 273 274 A = num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)]) 275 276 277 q_1 = operator.elliptic_multiply(q_in, include_boundary=False) 278 279 280 assert num.allclose( [-12.0-24.0/sqrt(5)], q_1.centroid_values ) 281 282 155 283 156 284 def test_mul(self): 157 operator1 = self.operator1() 158 operator1.apply_stage_heights(num.array([[1.0]])) #h=1 159 operator1.set_qty_considered(1) 160 operator1.build_boundary_vector() 285 operator = self.operator1() 286 287 q = Quantity(operator.domain) 288 q.set_values(2.0) 289 #q boundary_values should equal 0.0 290 291 operator.build_elliptic_boundary_term() 161 292 A = num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)]) 162 293 V1 = num.array([2.0]) #(uh)=2 163 294 U1 = num.array([[2.0],[0.0],[0.0],[0.0]]) 164 assert num.allclose(operator 1 * V1, 2*num.array(num.mat(A)*num.mat(U1)).reshape(1,))295 assert num.allclose(operator * q, 2*num.array(num.mat(A)*num.mat(U1)).reshape(1,)) 165 296 166 297 def test_cg_solve(self):
Note: See TracChangeset
for help on using the changeset viewer.